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ABSTRACT 

This work develops extensions and applications of a second-order land- 
surface parameterization, proposed by Andreou and Eagleson [1982], Procedures 
for evaluating the near surface storage depth to be used in one-cell land- 
surface parameterizations are suggested and tested by using the model. 

Sensitivity analysis to the key soil parameters is performed. A case 
study involving comparison with an "exact” numerical model and another sim- 
plified parameterization, under very dry climatic conditions and for two 
different soil-types, is also incorporated. 


i • 


1 : 
I ' 



ACKNOWLEDGEMENTS 


This work was completed with the support of the National Aeronautics 
and Space Administration (NASA) under Grant NAG 5-134. 

The work was performed by Stefanos A. Andreou, research assistant 
in Civil Engineering at llIT and was supervised by Peter S. Eagleson, Pro- 
fessor of Civil Engineering, 

Thanks are due to Antoinette DiRenzo, who performed all the neces- 


sary typing. 


TABLE OF CONTENTS 


PAGE 

TITLE PAGE 1 

ABSTRACT 2 

ACKNOWLEDGEMENTS 3 

TABLE OF CONTENTS ^ 

LIST OF FIGURES 6 

LIST OF TABLES 9 

NOTATION 10 

CHAPTER 1 Introduction 16 

CHAPTER 2 Sensitivity of Cumulative Yield During the Rainy 

Season with Respect to Changes of the Optimal Soil 
Properties K(l) and c 12 

CHAPTER 3 Sensitivity of Annual Yield to Changes of Storage 

Depth ■ 23 

CHAPTER 4 Sensitivity of Cumulative Yield to Storage Depth 28 

CHAPTER 5 Selection of the Appropriate Value of Storage Depth 36 

CHAPTER 6 Comparison with an Exact Numerical Model and Other 

Simplified Parameterizations Under Very Dry Condi- 
tions 38 

6.1 Introduction 38 

6.2 The Periodic Atmospheric Forcjiig 38 

6.3 The Parameterization of Fivicc-s in the Surface 

Boundary Layer 62 

6.4 The Soil Moisture Model 64 

6.5 The Force-Restore Method for Soil Temperature 

Prediction 35 

6.6 The Coupling with the Soil Moisture Model 68 

6.7 Evaluation of P^esults 20 


4 


PAGE 


CHAPTER 7 

Sunmiary and Conclusions 

80 


7.1 

Summary 

80 


7.2 

Conclusions 

81 


7.3 

Suggestions for Further Research 

83 

REFERENCES 



84 

APPENDIX 



86 


1. 

Program Taylor . Fortran 

87 


2. 

Program Winslow. Fortran 

101 


5 


LIST OF FIGURES 


Figure No, Figure Title PAGE 


1 Expected value of simulated annual yield as a func- 
tion of storage depth, Clinton, Massachusetts and 

Santa Paula California 26 

2 Comparison between simulated and observed CDF of 

annual yield, = 40cm, Clinton, Massachusetts 38 

3 Comparison between simulated and observed CDF of 

annual yield, = 100cm, Clinton, Massachusetts 39 

4 Comparison between simulated and observed CDF of 

annual yield, Z^ * 140cm, Clinton, Massachusetts 40 

5 Comparison between simulated and observed CDF of 

annual yield, Z^ » 160cm, Clinton, Massachusetts 41 

6 Comparison between simulated and observed CDF of 

annual yield, Z = 200cm, Clinton, Massachusetts 42 • 

r 

7 Comparison of storage change obtained from Milly’s 
and Eagleson's [1980] numerical model, Clinton, 

Massachusetts 

8 Comparison of total yield obtained from Hilly 's and 

Eagleson's [1980] numerical model, Clinton, Massa- 
chusetts 45 

9 Comparison between simiilated and observed CDF of 

annual yield, Z^ = 40cm, Santa Paula, California 46 

10 Comparison between simulated and observed CDF of 

annual yield, Z^= 100cm, Santa Paula, California 47 

11 Comparison between simulated and observed CDF of 

annual yield, Z^ = 140cm, Santa Paula, California 48 


Figure No. 
12 

13 

14 

15 

16 

17 

18 

19 

20 . 


Figure Title 


Comparison between simulated and observed CDF of 
annvaal yield, = 160cm, Santa Paula, California 

Comparison between simulated and observed CDF of 
annual yield, = 180cm, Santa Paula, California 

Comparison between simulated and observed CDF of 
annual yield, Z^ » 160cm, 30 years of simulation, 

Santa Paula, Califoraia 

Comparison of storage change obtained from Milly's 
and Eagleson's [1980] numerical model, Santa Paula, 
California 

Comparison of total yield obtained from Milly's and 
Eagleson's [1980] numerical model, Santa Paula, 
California 

Effective thermal conductivity as a function of 0 
and T. left; SILT LOAM. Right: Sand, After Milly 

and Eagleson [1980] 

Comparison of latent heat fluxes obtained by the 
numerical model (Milly and Eagleson, 1982) and those 
obtained by the analytical and numerical models when 
vapor flow is ignored. 

Liquid hydraulic conductivity, K and vapor conduc- 
tivity, D. as functions of 0 and T. Left: SILT LOAM 

Right: SAND, After Milly and Eagleson [1982] 

Moisture Retention in the SILT LOAM.Left; Main Branches 
of the hysteretic retention function. Right; Depend- 
ence of the first desorption curve on temperature 
(After Milly and Eagleson, 1982). 


PAGE 


49 


50 


51 


53 


54 


69 


73 


75 


77 


7 


Figure No. 


Figure Title 


PAGE 


21 


Comparison of latent heat fluxes computed from 
Milly's and Eagleson's [1982] mu-ierical model, 
Milly's and Eagleson's [1982] simplified one-cell 
parameterization and the analytical model, SILT 
LOAM 


78 


22 


Comparison of latent heat fluxes computed from 
Milly's and Eagleson's [1982] numerical model, 
Milly's and Eagleson's [1982] simplified one-cell 
parameterization and the analytical model, SAND 


79 


8 


LIST OF TABLES 


Table 

No. 

Table Title 

PAGE 

TABLE 

2.1 

Climatic and Soil Bropel'ties of Clinton 
Massachusetts and Sant4, Paula, California 

18 

TABLE 

2.2 

Sensitivity of Yield Due to Changes in Soil 
Properties and (Clinton, Massachusetts) 

20 

TABLE 

2.3 

Sensitivity of Yield Due to Changes in Soil 
Properties and (Santa Paula, California) 

21 

TABLE 

3.1 

Statistical Properties of Storm Characteristics 

24 

TABLE 

5.1 

Climatic and Soil Parameters to Calculate the 
Penetration Depth. 

55 

TABLE 

6.1 

Representative Values of the Forcing Parameters 
for Winslow in July 

60 

TABLE 

6.2 

Volumetric Heat Capacity of Soil Constituents 
(After DeVries, 1966) 

68 

TABLE 

6.3 

Climatic and Soil Parameters Winslow, Arizona 

71 


9 


NOTATION 


rJBSK TC 3t2 . i'nr Tj rt 


Symbols 

A 

c 


Definition and Dimensions 


d 

D 

e 

D 


ipv 




albedo 


[-] 


-3 -1 

volumetric heat capacity of a homogeneous medium [cal L deg ] 


coefficient for the sensible heat transfer 


coefficient for the water vapor transfer 


[-] 


[-] 


2 -2 -1 

specific heat of water vapor at constant pressure [L T deg ] 


pore disconnectedness index 


dlffusivlty index 


diaorption diffuslvity 


Vapor conductivity 


annual potential evapo transpiration 


annual actual evapo transpiration 


[-3 


[-3 

[l^t"^3 

[lt"^3 


[L3 


^[L] 


£ 

E 


exfiltration parameter [-3 

evaporation rate 

average annual actual evapotranspiration rate 
average annual potential evapotranspiration rate [LT 
actual evapotranspiration rate 
potential evapotranspiration rate 
transpiration rate from vegetation 
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[lt"^3 

[lt’^J 


,-ii 


[lt'^3 

[lt“^3 

[lt“^3 


Symbols 


Definition and Dimensions 



J 

K(l) 

k(l) 

k 

k 

V 

k 

s 

L 

L 

M 



ni_ 


exfiltraticsn capacity of soil 
infiltration capacity of soil 
heat flux into the soil 
surface retention capacity 
sensible heat 

incoming shortwave radiation 
down long-wave radiation 

rainfall rate 

evapotranspiration efficiency 
saturated hydraulic conductivity 
saturated intrinsic permeability 
Von Karman's Constant 
plant coefficient 
soil thermal diffusivity 
latent heat of vaporization 
Monin-Obukhov length 
vegetal canony density 
equilibrium vegetal canopy density 
rainy season length 
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Symbols 


Definition and Dimensions 


”i 


m 


m 

N 

n 


P 


F. 


Pa 

q* 


^Vb 



mean time between storms 

mean annual precipitation 

mean number of storms per year 
mean storm Intensity 
mean storm depth 

pore size distribution index of soil 
relative thickness of the atmosphere 
cloud cover 

effective porosity of the soil 
precipitation rate 
annual precipitation 
mean storm Intensity 
atmospheric pressure 

saturated atmospheric specific humidity 

specific humidity of the atmosphere at screen 
elevation 

bulk Richardson number 
annual groundwater runoff 

annual surface runoff 

gas constant 
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[-] 

[LT"^] 

[L] 

[-1 

[-] 
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[-1 

[LI 

VnT^] 

1 -] 

{-] 

[-] 

[L] 

[L] 

[LV^deg‘^J 
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t 
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t 

r 


U 

a 

w 



g 


s 


s 


- 1 /"? 

exfiltration "desorptivity" [LT 
infiltration ’’sorptivity'* [LT 
average soil moisture at the surface layer [-] 
average annual soil moisture at the surface layer [-] 


soil moisture concentration at time k [-] 

average annual atmospheric temperature [deg] 

air temperature at screen height [degj 

near surface soil temperature [deg] 

deep mean soil temperature [deg] 
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storm duration [T] 

time between storms [T] 

wind speed [LT ] 

upward capillary rise from the water table [LT ^] 

average annual yield [L] 

cumulative yield [L] 


percolation rate 

average annual percolation rate 

surface runoff rate 

average annual surface runoff rate 
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[lt“^] 

[lt“^] 

[LT~^] 
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y(s) 


Z 


r 
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Z 

o 

a 


0 
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surface layer thickness (storage depth) 
screen height 
surface roughness 
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Chapter 1 


Introduction 

This work consists of extensions and applications of the second-order 
Budyko-type parameterization of landsurface hydrology proposed by Andreou 
and Eagleson [1982]. It is based on a one-dimensional, short-term water 
balance model which can be coupled with a thermal balance model to obtain 
estimates of moisture and heat fluxes across the land surface. More speci- 
fically the objectives of this study were: 

1. Perform a sensitivity analysis of the annual yield estimated by 
the model with respect to the soil properties k(l) and c, about their 
''optimal" [Eagleson, 1982] values for contrasting climates. 

2 . Evaluate the sensitivity of the yield to the selection of the 
Storage depth used by the model. 

3. Establish an analytical procedure for making the above evalua- 
tion. 

4. Propose a way of selecting the storage depth independently from 
calibrations using detailed numerical models. 

5. Test the model under very dry climatic conditions and compare 
the results with those obtained by the parameterization or Milly and 
Eagleson [1982], 
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Chapter 2 


Sensitivity of Cumulative Yield During the Rainy Season with Respect 
to Changes of the Optimal Soil Properties lc(l) and c 

In .',11 tests of the short-term water-balance model, that appear 
in the Technical Report No. 280 [Andreou and Eagleson, 1982], for the 
catchments of Clinto*^^ 1-iassachusetts and Santa Paula, California, the 
values of soil Intrinsic permeability k(l) and pore disconnectedness 
index c were set equal to their optimal values, as they were derived 
by applying Eagleson *s [1982] ecological optimality hypotheses. In 
this section, a sensitivity analysis is performed, in order to find out 
the effect on the prediction of cumulative yield, of varying k(l) and c 
from those optimal values. 

It is important to know how robust the results of the model are with 
respect to changes in those values of soil properties, since in reality 
many uncertainties will be encountered about their true value. 

Thus, the short-term water-balance model, as described by Andreou 
and Eagleson [1982], was again applied at the two contrasting climates 
of Clinton, Massachusetts and Santa Paula, California. All climatic and 
soil parameters remained unchanged, as given in Table 2.1 except the 
values of k(l) and c that were varied one at a time. 

The model was run for a period equal to the rainy season length. 
First, the value of k(l) was Increased by 20 percent from its optimal 
value, everything else remaining constant, and the percentage change of 
cumulative yield from that corresponding to the optimal k(l) was calcu- 
lated. The same procedure was also followed by Increasing c by 20 per- 
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Table 2.1 


Climatic and Soil Properties of Clinton, 
Massachusetts and Santa Paula, California 


Clinton, Massachusetts Santa Paula, California 


M = 

0 

0.912 

M 

0 

0.424 

e * 

P 

0.150 cm/day 

fl> 1 

*o 

II 

0.274 cm/day 


3 days 


10.42 days 

r 

0.32 days 

m “ 

t 

r 

1.43 days 


365 days 


212 days 

w/e » 
P 

0 

w/e = 
P 

0 

% ' 

109 

% “ 

15.7 

% ' 

94 cm 

\ ' 

54 cm 

k “ 

V 

1 

k 

V 

1 

f 

a 

8.4^C 

f 

a 

13.8°C 

K = 

0.50 

K = 

0.25 

X 

0.578 

X 

0.0732 

k(l) = 

5.57xl0”^^*cm^ 

k(l) = 

12.27xl0“’^^cm' 

c » 

4.75 

c =* 

5.25 

n = 

0.35 

n = 

0.35 


[The values of M , k(l) , and c were set equal to those corresponding 
o 

to peak climatic values, according to the vegetal equilibrium hypothesis 
and the ecological optimality hypothesis, as they are described by Eagleson 
[1982.]] 
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cent from its optimal value and keeping everything else constant. It 
must be pointed gut that when k(l) or c was changed from its optimal 
value, the annual water balance [Eagleson, 1978] was solved in order to 
determine a new value of the annual average soil moisture s^ around 
which to linearize the evaporation and yield functions. 

The test was repeated using different values of (the surface 
layer thickness) in the range of 20cm 120cm. Two different tests 
were performed for each climate; one assuming bare soil and one by 
setting the vegetation equal to its optimum value M*. 

The results are shown in Tables 2.2 and 2.3. As was expected, for 
the humid climate of Clinton, Massachusetts the cumulative yield is 
rather insensitive to changes in the soil properties k(l) and c. There 
is no particular sensitivity to either of those two parameters. Also, 
there is very small difference between bare and vegetated soil. The ex- 
planation for this is that in the humid climate, evaporation is almost 
always at the potential rate (and for the linearized model used here, 
this happens all the time due to it’s structure). Thus, the only changes 
that can occur in the yield by varying the soil properties will be due 
to changes in storage. So, as k(l) increases less water is stored in 
the layer of depth Z^ and as c Increases more water is held in storage. 
The differences between bare and vegetated soil are very small and are 

due to small mnnerical differences between the functions of J(s , M, k ) 

o V 

used for bare and vegetated soil respectively. 

For the catchment of Santa Paula, we first observe a difference be- 
tween bare and vegetated soil. That is, in the presence of vegetation, 
control passes to the soil for longer time periods, so the role of evapo- 
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Table 2.2 


Sensitivity of Yield Due to Changes in 
Soil Properties and 



Clinton, Massachusetts 

(Bare Soil) 


s =0.71 
0 

s =0.75 
0 

z 

r 

(m) 

Z increase of yield 
due to 20% increase 
of U(l) 

% reduction of yield 
due to 20% Increase 
of c 

0.2 

11.11 

0.78 

0.4 

4.97 

2.72 

0.6 

5.27 

2.61 

0.8 

4.89 

3.58 

1.0 

4.10 

4.47 

1.2 

7.56 

2.49 


Z 

r 

(m) 

Clinton, Massachusetts 

s =0.71 
o 

% increase of yield 
due to 20% increase 
of k(l) 

i (M*«0.912) 

- o 

s =0.75 
0 

% reduction of yield 
due to 20% increase 
of c 

0.2 

3.17 

1.06 

0.4 

0.98 

2.43 

0.6 

2.28 

2.39 

0.8 

2.37 

3.42 

1.0 

1.85 

4.80 

1.2 

5.37 

2.34 


Table 2.3 


Sensitivity of Yield Due to Changes in 

Soil Pronerties and Z 
‘ r 

Santa Paula, California (Bare Soil) 


s *0.52 s =0.58 

o o 


z 

r 

(lo) 

% increase of yield 
due to 20% Increase 
of k(l) 

% reduction of yield 
due to 20% increase 
of c 

0.2 

5.07 

8.19 

0.4 

4.72 

20.33 

0.6 

8.76 

27.85 

0.8 

3.68 

23.36 

1.0 

2.79 

22.16 

1.2 

1.35 

21.45 


(m) 

Santa Paula, California 

L (M**0.424) 

s =0.55 
0 

% increase of yield 
due to 20% increase 
of k(l) 

s =0.60 
0 

% reduction of yield 
due to 20% increase 
of c 

0.2 

4.51 

9.72 

0.4 

12.16 

18.47 

0.6 

18.71 

31.45 

0.8 

11.95 

24.52 

1.0 

10.87 

22.70 

1.2 

11.69 

22.64 
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ration is reduced and percolation becomes more important than when M*0. 
Thus, the yield becomes more sensitive to changes in lc(l) in the presence 
of vegetation. 

We also observe that yield is more sensitive to changes in c than 
to changes in k(l) for both cases (M=0 and M*M*) . 

In both climates it turns out that knowing the true values of the 
soil properties k(l) and c is important in order to determine the soil 
moisture level in the layer near the surface. But for the humid climate, 
the accuracy of the estimates of those parameters does not significantly 
influence the estimates of the annual yield obtained by the model. On 
the contrary, for the semi-arid climate, deviations from the true values 
of k(l) and c can cause serious errors in the estimation of the annual 
yield . 


Chapter 3 

Sensitivity of Annual Yield to Changes of Storage Depth Z 

In Chapter 7 (Section 7,2) of Andreou and Eagleson [1982] a sensi- 
tivity analysis of the annual yield derived from the model was performed 
with respect to changes in the parameter Z^. The simulation period used 
was equal to the length of one rainy season (365 days for Clinton, tlassa- 
chusetts and 212 days for Santa Paula, California). The differences be- 
tween the statistics of the generated rainstorm events and interstorm 
periods and of the historical data v;ere shown in Table 6.1 of Andreou and 
E \gleson [1982]. Here the model was run for a longer simulation period 
corresponding to 10 consecutive years of successive precipitation events 
and dry periods. The statistics of the generated events and the corres- 
ponding historical values are shown in Table 3.1. A very small discrep- 
ancy between the two is observed. 

For every value of Z^ from 20cm to 200cm (using 20cm increments) the 
value of the average annual yield over the 10 year period was calculated. 

The precipitation characteristics that were used to generate the rainy 
and dry periods in Clinton, Massachusetts were those of Boston, Massachu- 
setts appropriately transformed, so that they corresponded to those of 
Clinton, Massachusetts. That was necessary to be done, since observations 
of annual yield, necessary for later comparison did not exist for Boston 
and on the other hand, hourly precipitation data from Clinton were not avail- 
able for analysis. 
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Table 3,1 


Statistical Properties of Storm Characteristics 


Storm depth 

Clinton, Massachusetts 
Historical 

(Boston, Massachusetts) 
E[h] » 0.86 

Generated 
E[h] = 

0.88 

[cm] 

Var[h] * 

1.50 

Var[h] » 

1.35 

Storm duration. 

E[t^] » 

0.32 

E[t I » 

0.32 

[days] 

r 

Var[t^J » 

0.10 

■ r 

Varft ] ■ 
r 

0.10 

Time between 

Elt^J - 

3 

ECt^J - 

3.11 

[days] 

VarCt^J - 

9 

Var[t^] ■ 

9.49 



Santa Paula, 

California 




Historical 

Generated 


Storm depth 

E[h] * 

3.41 

E[h] » 

3.31 

[cm] 

Var[h] = 

* 46.65 

Var[h] = 

37.35 

Storm duration 

E[t ] = 
r 

1.43 

E[t^] = 

1.49 

[days] 

Var[t ] 
r 

= 2.04 

Var[t^] “ 

2.38 

Time between 

E[t^] - 

10.42 

Eltbl - 

10.72 

[days] 

Var[t^l 

- 108.58 

Var[t^] - 

107.89 
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By using the sampled mean and variance of annual precipitation at 
Clinton, Massachusetts obtained from 30 years of data and assuming that 
the number of storms at Clinton is approximately the same as that in 
Boston (m^*l09) we can obtain the values of the parameters K and X of 
the rainstorm characteristics in Clinton by using the following formulas 
[Eagleson, 1978] 




[i + i] 


(3.1) 


and 


where 

nijj - mp /m^, nip * 111.3 
A A 

The derived values for K 
tively . 

Thus, using those values 
tation characteristics between Boston and Clinton to be the same, rainstorm 
events and interstorm durations were generated. 

It was found that for the humid climate of Clinton, Massachusetts, the 
value of remained always almost constant at 54cm, for any value of 
in the range 40 'U 200cm. On the contrary, for the semi-arid climate of 
Santa Paulas California, there is a drastic change of as varies, 
which can b^ seen from Figui*e 1. More precisely, there is a rapid de- 
crease of Y. as Z increases, although the percentage change of the yield 

A IT 

is reduced as becomes larger. 


"Hi - I 

2 2 

cm, Oy. » 268.38 cm , =» 109 

r , V 

A 

and X at Clinton are 0.73 and 0.175 respec- 
for K and X and assuming all other precipi- 
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EXPECTED VALUE OF SIMULATED ANNUAL FIELD Cctn3 



Expected value of simulated annual yield as a functicn of storage 
depth, Canton, Massachusetts and Santa Paula, California 
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i.' 0 

I 







\ *■ ' 


i ' 





The declining value of E[Y^] with increasing may be understood 
qualitatively by first realizing that the larger the "bucket" Z^,, the 
smaller will be the fluctuation in soil moisture within the bucket due 
to a given climatic forcing. 

In a dry climate where the evaporation rate is controlled by the 
soil, this reduced decay of soil moisture concentration during an evap- 
oration period means that the volume of evaporative flux increases over 
that for a smaller bucket having the same average soil moisture. This 
causes a redu‘ation in the waiter yield. 

Of course the average soil moisture concentration in the bucket is 
Itself dependent on the bucket size which may upset the above reasoning 
in a particular case. 

In a wet climate, the evaporation rate is climate controlled and this 
sensitivity is not present. 

It thus becomes important in dry climates at least to correctly de- 
fine Z . 
r 

In Chapter 4 a quantitative analysis is performed in order to explain 
the functional relation between the yield and the parameter Z^. Analy- 
tical expressions, relating the cumulative yield to Z^ are derived and the 
trends and behavior of they appear in Figure 1, are explained 

by using approximate solutions of those expressions. 


27 


Chapter 4 


Sensitivity of Cumulative Yield to Storage Depth 

In this chapter an attempt is made to develop an analytical rela- 
tion between the cumulative yield and the parameter during a precipi- 
tation and an evaporation event. The impact that variations of Z^ have 
on the yield is investigated by examining the sign and the order of mag- 
nitude of the derivative of the yield with respect to Z^. 

Following the development described in Chapter 4 of Andreou and 
Eagleson [1982], the one-dimensional short-term water balance of a soil 
column of depth can be written in the form: 

aZr ^ - i - 6 t - y (4.1) 

where 

e^ =» evapotranspiration rate 

y * yield rate 

i * rainfall rate 

n » effective porosity of soil 

“ storage depth 

s * average soil moisture concentration within the layer of 
thickness Z 

r 

By linearizing the values of e^ and y around their annual average values. 
Equation (4.1) can be written in the form: 

(i) During a precipitation event (assuming e^ ~ 0) 

nZ ^ = i “ Y C(s-s ) (4.2) 

r dt o 

(il) During an evaporation event (i*0) 
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(4.3) 


nZ 


ds 
r dt 



C* (s-s 

o 


) 


where y is the average annual yield rate that equals the average surface 

runoff rate y plus the average percolation rate y , and C and C* are 
s g 

the linearization coefficients which are given by; 


C 


m^tnt 2 3 


C 


e + 
P 



Vt 


(4.4) 


(4.5) 


and 



mean annual precipitation 
mean number of storms per year 
mean storm duration 

average annual potential evaporation rate 

where J is the evapotranspiration efficiency 


where 


C 


2 


3(^s/p) 

3s 


s=s 

o 



ilZs/p> 

3s 


s=s 

o 


(4.6) 


(4.7) 


p = m_ /(m^-nij. ) 

A r 

s^ = average annual soil moisture concentration 
Expressions for calculating C^, C 2 , and are given by Andreou and 
Eagleson (1982]. 


29 


An analytical solution for the differential Equations (4.2) and (4.3) is 
now derived. We again distinguish between the following cases: 


Precipitation 

The solution of Equation (4.2) is given in the following form: 


s(t) 


i - y + C s 


J L 




nZ 


1 - e 


nZ 


+ s^ e 


(4.8) 


where 

t^ * time that precipitation starts 

s. = initial soil moisture at time t 
i o 

By using its linearized form, the yield rate y during precipitation can 

be written as follows: 


y(s) » y(s^) + C(s-s^) 


(4.9) 


The cumulative yield y produced from time t to time t can thus be 

c o 

P 

written: 


rt 


[y(s ) - G s + C.s(t)]dt 
o o 


By setting t^ = 0 we obtain: 


^ - y(s ) + C s )t 

Cp o . o o o 


(i - y(s^) + C s^) - s^*C 


- 


-Ctn 

nZ 

nZ 

nZ 

r 

r 

r e 


C 

C 


- 




(4.10) 
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By differentiating y with respect to Z , we obtain a quantitative 
measure of the change that will occur in the cumulative yield during a 
precipitation event if we vary Z^. Thus, by using Equation (4.10) we 
have: 


E = 

dZ 

r 



-Ct 

-Ct n 
nZ 

te ^ 



E ! 

f nz 1 

1 r 

1 - e 

1 J 


i - y(s ) + C s - s.C 
o o 1 

c 1 

Z 

r 



(4.11) 


It is evident from Equation (4.11) that depending on the relative mag- 

Gp 

nitude of the components appearing in it, the sign can be either 

r 

positive or negative. Here, an attempt will be made to evaluate this 
derivative for a particular value of Z^ by assuming an average storm in- 
tensity and duration for the catchments of Clinton, Massachusetts and 
Santa Paula, Calfiomia. The chosen value of was 100cm. By using 
the soil and climatic properties of those two catchments given in Table 
2.1 we obtain: 


Clinton, Massachusetts 


i - y(s ) + C s - s. C = 2.68 -*0.5 + (7.69x0.72) - 7.69s. 
o o X X 

= 7.71 - 7.69 s. > 0, since s. < 1 
1 i — 

and 


-Ct 


f ^ 

nZ 

r -7.69x0.32'! 

te ^ 0.35 

, 0.35x100 

ll - e 

(. J 

Z 7.69 

r 

1 - e 


0.32 e 
100 


-7.69x0.32 

0.35x100 


Substituting in Equation (4.7), we obtain: 


SSStlSEl! 


Ti.- 
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dye 

-r^- -0.0001 [7.71 - 7.69 sj < 0 
r 


Santa Paula, California 


i - y(s ) + C s - Cs. * 2.69 - 0.35 + (4.32x0.55) - 4.32s. 
o o i i 


and 


4.71 - 4.32s, > 0, Since s. < 1 
X i — 


ri 

c 


1 - e 


nZ 

1 

J 


dv 

c 


-Ct 

nZ 


te 


-4.32x1.43 


0.35 
* 4.32 

= 0.0011 


1 - e 


0.35x100 


1.43 e 
100 


-4.32x1.43 

0.35x100 


Thus = -0.0011 [4.71 - 4.32sJ < 0 
aZ 1 


We observe that in both cases, at least on the average and in the 

vicinity of Z^ = 100cm, cumulative yield is expected to decrease as 

Z Increases, 
r 


Evaporation 

The solution of Equation (4.3) is given in the following form; 


s(t) 


1 




r "^‘^1 

('-e_ - y + C’s 1 


nZ 


fl - y + C s ’I 

nZ 

. T o 


1 - e ^ 

J, 

^ o 

1 - e ^ 

L C J 


X *“ e 

'T 

[ c 



C’(t -t 
0 


-Ct c'(t -t) 
o o 


nZ nZ nZ 

+ s, e *e ^ 

X 


(4.12) 
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where evaporation starts at time t^ ® t^, and the initial value of soil 

moisture at that time was calculated from Equation (4.8), by setting 

t » t . 
r 

The linearized function of the yield rate y during evaporation is 
given by: 

^g “ ^g^®o^ “ ®o^ 

Thus the cumulative yield during an evaporation event will be given by: 
t 

y. * C,.p.s(t) - C-.p.s Jdt (4.14) 

o r 


By substituting Equation (4.12) in Equation (4.14) » we obtain the 

following expression for y : 

*^E 


+ '=3” u 


+ C's^)(t-t^j' 

J 


c 


nZ 


nZ 


1 - e 


C 


-e„ -y +C's i-y+Cs 

7 . . g q 4. - Q 



r Gt 

o 

0 

ft 

o 

1 

rt 

1 

o 

rt 

o 


' 


nZ 

r 

ixZ 

r 

nZ 

r 


i 

“ y + C s 

o 

C' 

L, 

^6 

> 


. 

C ®ij 


(4.15) 


If we calculate the derivative of y^ with respect to Z^, we find 

E 

the following expression; 
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dy 


G. 


C^pn 


dZ 


nZ ’ C'(t -t) e 
, r . o 

^ ® 5T" 


O' (£„-th 


nZ 


C3pn 


- y + C 

's i - 


^ 

e 


-Ct 


r 0 


nZ 

r Ct 

e 

1 + 


1 


"Jol 


C j 

+ e 


nZ 


C'(t -t) ~ Ct 
o o 


nZ 


- 1 


ri - y + C S 


- s. 


(4.16) 


As we did for the precipitation event, here also the derivative of 
with respect to Z^ will be evaluated at Z^ * 100cm and assuming 


t « m. , t “ t 
r t^ o 


m. 


By substituting the parameters for Clinton, Massachusetts and Santa 
Paula, California in Equation (4.16), w:e get the following results: 


Clinton , Massachusetts 

cm 

We have: t^ = 0.32 days, t - t^ = 3 days, i » m^ * 2.69 C * 7.69, 

C' * 4.95, Z =* 100 cm 
* r 

dy 

= -0.040 - 0.030 (1.005 - s^) < 0, since < 1 
r 
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Santa Paula, California 


We have: * 1.43 days, t - t^ =* 10.42 days, i » = 2.48 C 

i"’-** 3.86, C^p “3.41, " 100cm 


We obtain: 



-0.273 - 0.114 (1.025 - s^) < 0, since < 1. 


4.32, 


Thus, in all cases we observe that at least on the average and for 
in the vicinity of 100cm, the cumulative yield decreases as in- 
creases. We also obser';e that at Santa Paula, California the yield 
is much more sensitive to changes in Z than it is at Clinton, Massachu- 
setts. At Clinton, Massachusetts, the value of is very close to 
zero. 

Those analytically derived results are consistent with the results 
obtained by the model and shown in Figure 1. 
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Chapter 5 


Selection of the Appropriate Value of Storage Depth 

Up to this point, the model was tested by using an apriori selected 
value of the storage depth and different results were obtained by vary- 
ing its value. In this chapter a way of determining through compari- 
sons with the observed values of the annual yield is discussed. 

To a first-order approximation, the observed expected value of the 
annual yield can be compared with the expected value of the annual yield 

obtained by the model after operating it for a certain number of Simula- 

♦ 

tion periods. The value of Z^ could then be fitted, so that the two ex- 
pected yields match. 

This comparison was made for the catchments of Clinton, Massachusetts 

and Santa Paula, California and the results can be summarized as follows. 

In Clinton, Massachusetts the expected value of annual yield obtained from 

30 years of observations (1904 'v 1933) is * 55.4cm. The expected 

value of the annual yield after a 10 year simulation period was found to 

be: E[Y. 1 = 54. 30cm. That value was found to be almost exactly the same 

A 

for a range of values of Z^ between 40cm and 200cm. This result indicates 

the insensitivity of the expected annual yield to the value of Z^ for the 

humid climate of Clinton. This can be explained by the prevailing climate 

control conditions in this area, as it was argued in detail in Chapter 2. 

Nevertheless, the result does not help us to determine the appropriate value 

of Z for this catchment. 

' r 

For Santa Paula, California, the expected value of the observed annual 
yield is = 17.4cm, and the value of the expected annual precipita- 

tion is m_ * 54cm. The simulated value of the expected annual yield 
obs 
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obtained by the model, as we observe In Figure 1, is decreasing as 

increases and reaches a value of 19cm at * 180cm. The average annual 

precipitation produced by the generated rai. ^*^rm events is equal to 

Hip » 63.58cm which is considerably larger than the observed value 
gen 

for Santa Paula. This is due to the fact that the average number of 
storms per year is very small (m^«15.7) and also the variability of the 
rainy season length was not taken into account in those tests (i.e., T 
was set equal to its average value m^). 

Thus again, by only using expected values of the yield, we cannot 
obtain conclusive results for the appropriate value of Z^. 

A more accurate way of fitting Z^ to observations of yield would be, 
by comparing the observed GDF of the annual yield to the CDF of the annual 
yield obtained by the model through simulation. This type of comparison 
was performed for the catchments of Clinton, Massachusetts and Santa Paula 
California. It was considered that the best fit between the observed and 
simulated CDF was achieved if they had similar shapes and slopes. Possible 
over or under estimations of the yield by the simulated CDF are expected 
due to the finite length of the simulation. 

In Figures 2 to 6 the values of the observed and simulated CDF's of 
the annual yield at Clinton, Massachusetts are plotted, for values of Z^ 
equal to 40, 100, 140, 160, 200cm, respectively. The precipitation char- 
acteristics of the rainstorm events were those at, Clinton using the derived 
values of k and A from Equations (5.1) and (52), We can argue that the 
best fitting between the two is achieved at Z^ - 160cm where a very good 
agreement with the observed values of yield exists. This result strongly 
indicates that a value of Z^ in the vicinity of 160cm will be appropriate 
for operating the model. 
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By using * 160cm, a comparison between the fluxes obtained by 
the analytical model and the isothermal version of C. Milly's [1980] 
numerical model was performed. This type of comparison is described 
with more details in Andreou and Eagleson [1982], where = 100cm 
was assumed. The results for the storage change and for the yield 
produced, are shotm in Figures 7 and 8 respectively for Clinton, Massa- 
chusetts. It is evident that by using = 160cm, we clearly have an 
improvement of the analytical model. 

The same test was also performed for Santa Paula, California. The 
simulated and observed CDF's of the annual yield are shown in Figures 9 
through 13 for values of Z^ equal to 40, 100, 140, 160, 180cm, respectively. 
It appears that when Z^ is again about 160cm we obtain the best fitting 
between observed and simulated values of the yield CDF. Since the num- 
ber of storms per year in Santa Paula is small, we expect to obtain even 
better results if we run the model for a longer simulation period, since 
we will approach even closer the historical statistics of the precipita- 
tion events. We must nevertheless, keep in mind that some of the dis- 
civepancies between observed and simulated values of the yield are due to 
the. fact that the variability of the rainy season length was neglected 
during those simulations. 

For Z^ = 160cm, where the best fitting was observed, the model was 
operated for a longer simulation period equal to 30 years. The obtained 
CDF of the simulated annual yield is compared with the observed in Figure 
14. It can be argued that it gives a fairly good estimate of the actual 
CDF. 
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CLINTON, MA 

□ NUMERICAL MODEL 
O analytical model (2r=IOOcm) 
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FIGURE 9 

Comparison between simulated and observed CDF of 
annual yield, Santa Paula, California 
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figure 11 

Comparison between simulated and observed CDF of 
annual yield, Zj. = 140cm, Santa Paula, California 
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By setting Z^* 160cm, the fluxes obtained by the analytical model 
and by Milly’s [1980] isothermal version of his numerical model are now 
compared again for Santa Paula, California. Storage changes are shown 
in Figure 15 and yield produced is shown in Figure 16. It is apparent 
that setting » 160cm produces an improvement over the previously obtained 
results where » 100cm. 

The "best" value for Z^ obtained through simulation can now be com- 
pared with the value of the penetration depth of a step change in siirface 
soil moisture corresponding to the average climatic and soil conditions of 
the region under investigation. 

The value of the penetration depth, combining the diffusive component 
with the gravitational seepage component is given by [Eagleson, 1978]: 


Z 


max 


= 4(Dt)“ + 


tK(0^) 

n 


(5.1) 


where D is the sorption diffusivity (D^) or the desorption diffusivity 

(D ) , t = m if D » D . and t » m if D * D , K(© ) is the hydraulic conduc- 
6 tl X C« 6 O 

r b 

tlvity at the average soil moisture level 0^ and n is the effective poros- 
ity. 

Values of D, t, K(0) and n for Clinton, Massachusetts and Santa Paula, 
California are given in Table 5.1. By substituting in Equation (5.1) we 
find: 

Clinton, Massachusetts 

For infiltration: Z^ * 65.53 +0.33 = 65.86cm. 

For exfiltration: Z * 50.04 + 3.11 * 53.11cm 

e 
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FIGURE 15 

Comparison of storage change obtained from Milly's 
gleson’s [1980] numerical model. Santa Paula 
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FIGURE 16 

Comparison of total yield obtained from Mill’y 
ai^d Eagleson's [1980] numerical model, Santa Paula, California 

54 





Table 5.1 


ClinCon, Massachusetts 


^3 2 

m » 0.32 days , D. » 9.71x10 cm /sec 


-4 2 

m * 3 days , D = 6.04x10 cm /sec 

b ® 


K(l) - 4.20x10"^ — • " “ 0*35 


Santa Paula, California 


-2 2 

» 1.43 days , = 1.18x10 cm /sec 


m. » 10.42 days , D * 3,456x10 cm /sec 
t, ’ e 

b 


K(l) = 9.25x10"^ n = 0.35 
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Santa Paula, California 


Tor infiltration; Z. - 153 i- 3.26 * 156.26cm 

1 

For exfiltration: Z = 70.55 + 23 = 93.55cm. 

e 

Although the calculated value of the penetration depth almost coin- 
cides with the value obtained through simulation for the catchment of 
Santa Paula, there is a discrepancy of about Im between the two for the 
catchment of Clinton, Massachusetts. 

Observing Figures 2 to 6, we see that the selected value of at 
every simulation run, influences the fitting of the CDF's primarily dur- 
ing the very dry years at Clinton- Since the analytical model uses a lin- 
earization around the average annual soil moisture in order to estimate 
the evaporation flux, the tangent to the evapotranspiration efficiency 
curve at the point corresponding to has a very rMall slope. Thus, the 
values of evaporation predicted by the model when the soil moisture level 
becomes very low are considerably overestimated. This results in less 
water stored in the bucket oi depth Z^ and thus eventually in less percola- 
tion to the water table. This eTsplains the fact that by choosing a value 
of Z^ on the order of 60cm as predicted by the penetration depth calcula- 
tions an underestimation of the yield is obtained during the very dry years. 

On the other hand, the value of Z^ does not play a significant role 
during the wet years for a humid climate, because the evaporation rate is 
usually equal to the potential evaporation rate. In any case, the dif- 
ferences in the CDF's for the different values of Z^ are not that pro- 
nounced as in the semi-arid climate of Santa Paula- From the above obser- 
vations it is found that for a semi-arid climate, such as that of Santa 

Paula, California, the value of Z^ obtained through penetration depth 
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calculations is very close to that obtained by fitting the observed CDF 
of the annual yield to that predicted l>y the analytical model through sim- 
ulation. It is also found that the values of Z obtained by the above 
mentioned procedures differ significantly for the humid climate of Clinton. 
Fortunately, however, knowing the appropriate value of for a humid 
climate becomes c? some importance only for the case of very dry years. 

As it will also be shown in Chapter 6, the value of is highly vary- 
ing depending on both the climate and the soil of a region. Thus, it is 
not appropriate to select it arbitrarily as is done in some algorithms. 

For example, Budyko [1956] suggested a value of Z^ - 100cm, Arakawa [1972] 

assumed Z^ - 30cm, Gates et al. [1977] suggestes 0^ Z^ ~ 30cm, where 0^ 

^c c 

is the field capacity and Shukla [1977] proposed 0^ ~ 10cm. 

In summary, this research demonstrates that the important climate 
and soil conditions of a region can be incorporated into a priori estima- 
tion of Z^ through computation. of the penetration depth. T-Hiere long-term 
water yield data are available Z^ may be estimated by fitting simulated to 
observed cdf's of annual yield. 


Chapter 6 


Conparlson with an Exact Munerical Ilodel and Other Sinplifled 
Parameterizations Under Very Dry Conditions 

6.1 Introduction 

The latent heat flux obtained by the model is compared with that 
obtained by the numerical model for heat transfer and moisture flow in 
soils developed by Milly and Eagleson [1980]. A comparison is also made 
with other simplified parameterizatlons proposed by Milly and Eagleson 
[1982]. 

The climate chosen to demonstrate this comparison is that of Winslow, 
Arizona, which is characterized by very dry conditions. The model was 
tested for two types of soil; silty loam and sand. A periodic atmospheric 
forcing was applied for a period of ten days. The force-restore method 
was used in order to update the estimates of the near surface and the 
deep soil temperatures. Thus, a thermal balance model was operated con- 
junctively with the soil moisture model, the coupling between them occur- 
ing through the evaporation rate e^ and also through changes in the soi"'. 
emisslvity e and surface albedo A due to changes in soil moisture. 

6.2 The Periodic Atmospheric Forcing 

The surface boundary layer is forced by six atmospheric variables, 
which ares the incoming shortwave radiation 1^, the down long-wave radi- 
ation from the clouds I. , the precipitation rate P, the air temperature 

a 

T , the wind speed U , and the vapor pressure of the air p . 
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The type of periodic forcing chosen is the one described by Milly 
and Eagleson [1982] and it will be briefly repeated here for convenience. 

a. Shortwave Radiation 

The intensity of shortwave radiation reaching the ground is 
given by: 


I 

s 


fw, 


BO 


-n'a m 2 

sin a e (1 - 0.65N ) 


sin a > 0 
sin a < 0 


( 6 . 1 ) 


-2 -1 

where W-no the solar constant (2 cal cm min ) , a is the angle 
i>U 

of the sun above the horizon, a^ is the molecular scattering coef- 
ficient, n* is a turbidity factor, m is the relative thickness of 
the atmosphere and N is the proportion of the sky covered by clouds. 
The angle a is given by: 


sina » sin6 


sin({) + cos6*coscf) cos 



( 6 . 2 ) 


where <S is the solar declination, <f> is the latitude and t is the 
time in hours since midnight. 


a^ - 0.128 - 0.054 log m (6.3) 

and 

m = (sina)'*^ (6.4) 

Representative values of the forcing parameters for Winslow in 
July, which were used in all applications of the model to be described 
later, are shown in Table 6.1. 
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Table 6,1 


Representative Values of 

the Forcing parameters 

for Winslow 

in July 

Parameter 

Winslow/ July 

^a 

12.9*C 

P 

18.6 cm yr ^ 


35® 01' 

6' 

21® 30' 

n’ 

2.5 

N 

o 

0.45 


94 hours 

t 

r 

2 hours 


20 hours 

X 

0.193 cm hr“^ 

T 

m 

25.8®C 

"di 

7.8®C 


15 hours 

U 

m 

360 cm s ^ 

“di 

180 cm 3 

t 

u 

18.5 hours 

Pva 

7 in-6 -3 

7x10 g cm 


60 


The cloud cover ratio N is taken as follows: 


N if p = 0 

N » ] ° 

1 if p > 0 

k 

b. Down Longwave Radiation 

The atmospheric longwave radiation is given by: 
Ig, » £ a(T + 273)^ (1 + 0.17N^) 

jCd d Si 

where the atmospheric emissivity is 

e = 9.37x10"^ (T + 273)^ 
a a 

c. Precipitation 


(6.5) 


( 6 . 6 ) 


(6.7) 


The precipitation rate P is expressed as a periodic function 
of time of the following forms: 


0. ti - + K(t^ + t^) < t < + K(t^ + 

1, t^ + K(t^ + tp < t < K(t^ + tp + t^ 


( 6 . 8 ) 


where i represents the average storm intensity, t and t are re- 

b r 

presentative values of the time between storms and storm duration 
respectively, t^ is the starting time of the first storm, and K is 
any integer. Values of the parameters appearing in Equation (6.8) 
are given in Table 6.1. 

d. Air Temperature and Wind Speed 


Monthly averages of the three-hourly, diurnally varying air 
temperature and windspeed for Winslow were fitted to the following 
cosine curves: 
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■=“= 'iff 

cos (t - t_^)] (6.10) 

X 


Values for T , T, , t_, U , U, , and t are given in Table 6.1. 
m d^’ T* m’ d^’ u ® 

The value of p was assumed constant and is given in Table 6.1. 
'^va ® 

6.3 The Parameterization of Fluxes in the Surface Boundary Layer 


a. Potential Evaporation Rate 


The evaporation rate is calculated through the aerodynamic 
Equation: 

c U 

E (p _ p ) 

P-. '^va vg 


( 6 . 11 ) 


where c^ is a bulk transfer coefficient, p^ is the density of liquid 
water, p^^ is the density of water vapor at the ground surface, and 
p^^ is the density of the water vapor in the air. 

Equation (6.11) was used in the model only to evaluate a change- 
ing value of the potential evaporation rate. When the surface became 
dry, the evaporation rate £ was calculated by the Equation: 


E » min(e , e ) 

T p 

where: e^ = e^ + 6^(3 - s^) 

which was documented by Andreou and Eagleson [1982], and where e^ is 
the annual average evaporation rate,.ep is the annualxpotential evap- 
oration rate and is a linearization coefficient. 


b. 


Sensible Heat Transfer 
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The sensible heat transfer was expressed by: 


H » -p c U (T - T ) 
p H a' a 


(6.13) 


where p is the air density, c is the specific heat of water vapor 

« p 

at constant pressure, c is a oulk transfer coefficient, T is the 

n a 

air temperature and T is the ground temperature. 

O 

Under conditions of neutrisl stability, the transfer coefficients 
become; 






Z >>2 
iln ^ 
o-' 


(6.14) 


where k is Von Karman’s constant (=0.4), Z is the screen height and 
is the surface roughness. 

Under unstable conditions the transfer coefficients are func- 
tionally related to their neutral values through [Anderson, 1976]. 



2 tan ^(x) 



(6.15) 


where 

X 

X = (1 - 16 Z /L) (6.16) 

a 

and L is the Monin-Obukhov length, which is related to the bulk 
Richardson number 
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2 Z (T - T ) 
g a 

(T + T ) U ^ 
a g a 
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(6.17) 


through the expression: 


'H 








1 - 




(6.18) 


Jin 


1+x 


+ Jin 


1+x 


2 tan”^(x) + ^ 


Equations (6,15) - (6.18) were solved iteratively and a table rela- 
ting (R.)„ to c„ and c was created. 

1 o rl W 

For stable conditions the following relation holds: 


'H 


w 






1 - 


»i>B 


• <Vb " *1 


cr 


cr 


(6.19) 


, (Rj) 3 i R^ 


cr 


where R. is the critical Richardson number, equal to 0.2. 
cr 

In all applications that follow » 200cm, 0.1cm and 


<=H>N ■ 

6.4 The Soil Moisture Model 

The model for updating the soil moisture within a surface layer of 
thickness Z^, is the one developed by Andreou and Eagleson [1982], and 
the linearized equations of the short-tem water balance were given in 


Chapter 4 (Equation 4.2 and 4.3). The only difference here is that dur- 
ing a precipitation event (Equation 4.2), the evaporation rate, is set 
equal to the changing potential rate e , in order to be consistent with 
the parameterization of Milly and Eagleson [1982] with which the com- 
parison Is made. 
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In order to locate more accurately the passage from climate control 
to soil control after a precipitation event, it is necessary to calculate 
the time t^ after precipitation ceases at x^hich the surface becomes dry. 
This time is given by: 



( 6 . 20 ) 


where the desorptivity is given by: 


S 

e 


2 s 


l+d/2rnK(l)i/^(l)c|)^(d)-i- 
ixm 


( 6 . 21 ) 


and s is the average soil moisture concentration within the layer of 
thickness 2^, imtutdlately after the precipitation ends. For a more de- 
tailed reasoning of this procedure see Andreou and Eagleson [Chapter 7, 
Section 7.5, 1982]. 

6.5 The Force-Restore Method for Soil Temperature Prediction 


The linear differential equation for estimating the surface tempera- 
ture Tj^ is given [Deardorff, 1978] by: 


^ G - - T^) (6. 

The values of c^^ and C 2 are given by: 



where A is the thermal conductivity and C is the volumetric heat capacit 3 r 
of the homogeneous medium. 
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The heat flux into the soil G, is expressed from a surface energy balance as 


G - (1 - A)!^ + E - O^CL + T^)E - H + Ojj I^P 


(6.25) 


where A is the surface albedo, the incoming shortwave radiation, £ 
the surface emmissivity, the downward longwave radiation, I 9 ,, the 


Hn 


upward longwave radiation (» ea(T + 273) ), E the evaporation rate, 

I z*o 

and c^ the density and specific heat of water, H the sensible heat 
loss, T^ the air temperature and P the precipitation. The heat: loss due 
to surface runoff and surface detention storage are neglected as not 
important. 

The deep soil temperature T 2 , which varies slowly due to the annual 
cycle of forcing is obtained from [Deardorff, 1978] 


fa 

dt 


(A Q N, T)“^*G 
a 


(6.26) 


The value of used in the simulations described by Milly and Eagleson 
[1982] was set equal to 20. For the reasoning behind this, see Milly and 
Eagleson [1982, Section 4.4]. 

6.6 The Coupling with the Soil Moisture Model 

The coupling between the thermal and water balance models occurs 
not only through the value of the evaporation rate E which was discussed 
in more detail earlier, but also through changes in the moisture content 
which influences the albedo, the emissivity, the thermal conductivity 
and the heat capacity of the soil. 
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Since the soil moisture model used does not predict the (volumetric) 
soil moisture at the surface 0^, but an average soil moisture 0 within 
the layer of thickness Z^, the value of 0^ will be approximated as 

fo if e_ < e and t > t 

0, - P ° (6.27) 

[5 If t < t 

where t is the time passed after precipitation has ended. 

The same t 3 rpe of approximation is also made by Milly and Eagleson 
[1982, Section 4.4.2] in their parameterization. 

The volumetric heat capacity of soil C is expressed as a weighted 
average of the capacities of its components [de Vries, 1966]: 

5 

C » I c. 0. (6.28) 

i*l ^ 

where 0, and c. are the volumetric fraction and the volumetric heat 
i 1 

capacity of the i'th soil constituent. The five soil components are (1) 
water, (2) air, (3) quartz, (4) minerals, (5) organic matter. The heat 
capacity of each constitutent is given in Table 6.2. The volumetric frac- 
tions for silty loam and sand were given in Table 6.3. The effective ther- 
mal conductivites X for silty loam and sand, as a function of 0 and T 
were calculated by Milly and Eagleson [1982] and are shown in Figure 17. 

The product AC appearing in Equation (6.23) and (6.26) of the force- 
restore method was evaluated in the manner described by Milly and Eagleson 
[1982] and will be repeated here for convenience, (The subscript *’2" is 
used when we refer to the prediction equation for T^) • Thus, we have: 

(AO 2 » A(0.) c(Op (2.29) 
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Table 6.2 


(After DeVries, 1966) 

Volumettlc Heat Capacity of Soil Constituents 


Constituent 

i 

!i 

liquid water 

1 

1.0 

air 

2 

3x10' 

quartz 

3 

0.46 

other minerals 

4 

0.46 

org. matter 

5 

0.6 



FIGURE 17 

Effective thermal conductivity as a function of 0 ?nd 
left; SILT LOAM. Right; Sand, After Milly and Eagleson [1980] 




where 0^ is the Initial moisture content, since due to the short duration 
of the simulations performed, the departure from initial conditions will 
be small. 

T'Jhen Ac is evaluated for use in Equation (6.23) however it is impor- 
tant to account for the time-varying surface moisture. Thus we have: 

(AC)j^^ » 0.3 [A(0^) c(0^)]^ + 0.7 (AC) 2^ (6.30) 


where the subscript "1" refers to the prediction equation for and the 
value of 0^ is given by Equation 6.27. Equation (6.30) is the one ap- 
plied in Milly's and Eagleson's [1982j parameterizations and consists of 
a slight sisjplification of a procedure proposed by Deardorff [1978]. 

The value of changing albedo is calculated as follows: 


A 


r 2w 

/d ■■ - ^d> -r ” 

A 20 > n 

w 


(6.31) 


Values of A^ and A^ are given in Table 6.3. For the soil emmissivity £, 
we will use a value equal to 0.95 if 0^^ 0 and a value of 0.9 if 0^ = 0. 

The s6il-moisture and the force-restore equations were solved simul- 
taneously using an explicit numerical procedure. The time step of inte- 
gration was equal to a quarter of an hour. 

'6.7 Evaluation of Results 

The latent heat flux calculated by the proposed parameterization 
was compared with t’at obtained by using the numerical model developed 
by Milly and Eagleson [1980] and with other simplified parameterizations 
proposed by Milly and Eagleson [1982] . The climatic variables and soil 

properties used are shown in Table 6.3. 
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Table 6.3 


Climatic and Soil Parameters 


Winslow, Arizona 


e = 

P 

0.449 cm/ day 


4.58 days 

®t " 

0 . 10 days 

r 


m^ = 

365 days 

h 

0 

0.1 cm 

w/e s* 
P 

0 

T 

a 

12.9“C 

“ 

74 

m ■ 

Pa 

22.33 cm 

K « 

0.32 


For Silty Loam 
n * 0.46 

k(1) s» 1. 24x 10”^ cm^ 
c 5 
- 0.20 

A = 0.10 
w 


For Sand 
n = 0.35 

k( 1) = 2.48x10 ^ cm^ 
c 5 

A, = 0.35 
a 

A « 0.25 
w 


I 

i 

t 

f 
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The simulation period lasted 10 days and the initial soil tempera- 
ture was set equal to 24°C while the initial soil moisture concentration 
was set equal to 0.083 for the silty loam and 0.091 for the sand. The 
results for the silty loam are shown in Figure 18. The solid line repre- 
sents the solution obtained using the ’’exact” numerical model. This model 
is fully documented by Milly and Eagleson [1980] . The open circles repre- 
sent the solution obtained using the numerical model, but where vapor flow 
is neglected as described by Hilly and Eagleson [1982; Section 6]. The 
black dots represent the results obtained using the parameterization de- 
scribed in this report. The value of the storage depth was fitted in 
order to obtain the best approximation with the numerical model. The op- 
timal value of Z^ was found to be equal to 2.97cm. We observe that the 
estimates of latent heat are in very good agreement with those of the 
numerical model up to the time that control passes to the soil. This 
occurs about nine hours after the end of the precipitation. We observe 
that when control passes to the soil there is a sudden drop in the latent 
heat flux, which is now considerably less than the one predicted by the 
numerical model. The reason for that is that vapor flow plays an impor- 
tant role in the early stages of exfiltration, when control passes to 
the soil and also when the soil moisture level in the surface layer is 
very low, as it is here. The effect of neglecting vapor flow in such a 
case can be seen very clearly from the solution of Hilly's [1982] numer- 
ical model, as it is plotted with the open circles in Figure 18. Here 
the plotted circles represent the value of latent heat by the numerical 
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model when vapor flow is neglected. We observe that the results are very 
similar to those of our parameterization (as shown by the solid circles 
in Figure 18) , where vapor flow is not included in the equations modeling 
the moisture dynamics. The importance of vapor flow under very dry condi- 
tions can also be seen from Figure 19, where we observe that the values of 
hydraulic conductivity K(0) and vapor conductivity are of the same 

order to magnitude when 0 is in the vicinity of 0.1, which is the case in 
our experiment. It is evident that vapor flow will be important only under 
very arid conditions and for particular types of soil, as is easily seen 
from Figure 19. If we want to take vapor flow into account in our model, 
a modification is necessary. As is proposed by Ililly and Eagleson [1982] , 
an effective value of the diffuslvity can be calculated, in which vapor 
flow is explicitly considered. The exfiltration capacity of the soil 
can then be evaluated through the selection of an appropriate formula. 

Here, we will evaluate the exfiltration capacity by using the Philip [1960] 
equation 


f (t) = j S t"“ - -5-[K(0.) + K(0 )] 
e ^ e z L o 

where 


2 ( 0 - 0 ) [- ^ 


D^ce.c) j 


and 


D (0, t) » 1.85 S 


-1.85 f* _ g^^^^jO.85 


K[0g(*)l + } 


dij; 


(6.32) 


(6.33) 


(6.34) 



Left: SILT SOAM Right: SAND. After Milly and Eagleson [1982] 


Equation (6.32) was applied to estimate the evaporation rate only 
from the time that the surface became dry (0j^»O) , up to the time when 
the evaporation rate obtained by (6.32) was of the same order of mag- 
nitude with the evaporation rate obtained from our original model. As 
we see from Equation (6.34) the effective diffusivity is updated at every 
time step. The integral appearing in Equation (6.34) is approximated 
using the functions of hydraulic conductivity, vapor conductivity and 
matrix potential, shown in Figures 19 and 20. 

The results obtained when this procedure is applied, are shown in 
Figure 21, by plotted circles. As we see those results are in almost 
perfect agreement with those obtained by Milly’s and Eagleson's [1982] 
parameterization. It should be noted that the computational burden in- 
troduced by these modification does not exceed that of Milly's simplified 
parameterization, although depends upon the form of the K(0) , 
and ip(0) functions chosen. 

The results obtained for the sandy soil are shown in Figure 22. 

Vapor flow was neglected in this case, since as we can see from Figure 
19, the vapor conductivity is much smaller than the h3;^draulic conductivity 
for the sandy soil and for 0 > 0.01. Again here we observe fairly good 
agreement between the proposed parameterization and that by Milly and 
Eagleson [1982] . The initial discrepancies of both from the numerical 
solution, are due to a transient error because of non-equivalence of 
initial conditions. In this case, the optimal value of the was found 
to be equal to 15 cm. 
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Moisture Content Moisture Content 

FIGURE 20 

Molstuxe Retention in the SILT LOAM. Left: Main Branches of the hysteretic 

retention function. Right: Dependence of the first desorption curve on 

temperature (After Mllly and Eagleson, 1982) 
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FIGURE 21 

Comparison of latent heat fluxes computed from Mllly's and Eagleson's 
11982] numerical model, Mllly's and Eagleson’s [1982] simplified one-cell 
parameterization and the analytical model, SILT LOAM 
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Chapter 7 


7 .1 Summary 

In this study, extensions and further applications of a second-order 
Budyko-type landsurface parameterization [Andreou and Eagleson, 1982] were 
made, 

A sensitivity analysis was first performed, in order to examine the 
changes of the cumulative yield predicted by the model, caused by changes 
of the soil properties k(l) and c from their values derived by applying 
the ecological optimality hypotheses. [Eagleson, 1982] 

A procedure was developed for obtaining analytically the sensitivity 
of the yield to the near-surface storage depth as defined by the model. 
Thus, for the proposed model and by using the soil and climatic proper- 
ties of a given region, it is possible to derive analytically a measure 
of the sensitivity of the yield to the near surface storage depth, 

A methodology of assessing the ’’best'* value of the storage depth is 

proposed. The CDF of the annual yield obtained through simulation by the 

model was compared to the CDF of the observed values of annual yield, and 

the value of storage depth that gave the best fitting between the two was 
selected. The validity of the above methodology was tested through compar- 
isons of the results obtained by the model and those obtained by applying 
Willy's [1980] numerical model. Having established a value of the storage 
depth by applying the previous method, the results of the comparison were 
always better than those obtained through setting the storage depth at its 
"nominal" value of Im (as is suggested by several Investigators) . 
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Finally, comparisons of the results obtained by the proposed model 
were made with results obtained by Milly and Eagleson [1982] using other 
simplified parameterizations . The importance of vapor flow under very 
dry conditions and for certain types of soils was investigated. Necessary 
modifications of the model, in order to handle those conditions, were sug~ 
gested and tested. 

7.2 Conclusions 

The conclusions derived from this research are the following 

1. The annual yield obtained by the model, is sensitive to the 
values of k(l) and c derived from the ecological optimality hypotheses 
[Eagleson, 1982]. But for the humid climates, the accuracy of the 
estimates of those parameters does hot; affect significantly the es- 
timates of the yield. On the contrary, for the semi-arid climates, 

it was found that the yield was very sensitive to those parameters. 

For the tested climates, the yield was also found to be much more sen- 
sitive to the value of the pore disconnectedness index c than to the 
value of the saturated intrinsic permeability k(l) . 

2. For the two contrasting climates of Clinton, Massachusetts and 
Santa Paula, California it was found by the model and also verified 
analytically, that the yield was much more sensitive to the selected 
value of storage depth for the semi-arid climate than for the humid 


climate. 


3. If observations of the annual yield are available, the appropri- 
ate value of the storage depth to be used by the model, can be asses- 
sed through comparisons of the CDF's of annual yield, observed and 
simulated. The validity of obtaining the value of storage depth by 
applying this technique was verified, when the model was operated in 
real time and comparisons were made with Milly's numerical model. The 
possibility of a priori selecting by setting it equal to the penetra' 
tion depth was also indicated by this study. 

It was found, from one application, that for a semi-arid climate 
the value of determined with the above method was very close to the 
value of the penetration depth, thus providing the possibility fivr 
a priori selection of for such a climate. Although for a humid cli- 
mate, the same result was not found, it was established that knownlng 
the accurate value of Z^ is of importance in humid climates only dur- 
ing the very dry years. 

4. It was found that the model with its present structure could not 
handle extremely dry situations for certain types of soils, where 
vapor flow is Important during exfiltration. However, if the vapor 
conductivity dependence upon soil moisture is known, then they can 
be Incorporated into the model, as suggested in Chapter 6. If this 
is done, it was found that the model can give very statisfactory re- 
sults, when calibrated with Milly's [1980] numerical model. These 
results where very close to those obtained by the simplified parame- 
terization of Milly and Eagleson [1982] which is calibrated to the 
numerical model through a fitted moisture redistribution parameter. 
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5. From this research It was found that a wide range of the appro- 
priate value of storage depth exists in reality, depending on the type 
of climate and soil for every region. Thus, for one-cell models of 
landsurface parameterization, we must be careful in the selection of 
the storage depth. Choosing it to be uniformly equal to Im, as is 
very often done, can yield large errors in the computed surface fluxes. 

7 . 3 Suggestions for Further Research 

In addition to further tests to verify the model, more extensive re- 
search is needed, to study the interrelation of storage depth, c Umate and 
soil for a variety of climates and soils. 
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PROGRAM TAYLOR. FORTRAN 
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C THIS PROGRAM GENERATES RAINSTORM EVENTS. STORM DURATIONS 
C AND INTERSTORM PERIODS WHICH PRESERVE THE HISTORICAL STATISTICS. 

C IT CALCULATES THE SOIL MOISTURE OVER A DEPTH CLOSE TO THE SURFACE 
C EVERY HALF HOUR .IT PLOTS THE EVAPOTRANSPIRATION 
C FUNCTION. THE SURFACE RUNOFF AND PERCOLATION FUNCTIONS. IT ALSO 
C PLOTS THE DAILY SOIL MOISTURE DURING THE RAINY SEASON LENGTH. 

C IT CALCULATES THE TOTAL STORAGE CHANGE, THE CUMULATIVE 
C EVAPORATION AND YIELD AT THE END OF EVERY RAINY OR 
C INTERSTORM PERIOD 

C IT HAS THE OPTION OF USING MANABE'S MODEL 
C TO CALCULATE THE MOISTURE FLUXES 

C THE VALUE OF THE POTENTIAL EVAPOTRANSPIRATION RATE 
C IS SET EQUAL TO ITS ANNUAL AVERAGE VALUE 
C THIS PROGRAM ALSO CALCULATES VALUES OF 
C ANNUAL YIELD FOR A SPECIFIED SIMULATION 
C PERIOD GREATER THAN ONE YEAR. 

C CLIMATIC AND SOIL VARIABLES 

c epp=averag9 annua] evapotranspiration ‘rate(cm/day) 
cmtb=mean t]me between storms (days) 
c njtr=mean storm duratlon(days) 
c mpa=mean annual predpl tatlon(cm) 
c mtau=mean rainy season length(days) 
c ta*average annual a<r temperature(C) 
c mnu=mean number of storms per year 
c n=soi1 porosity 

c k 1=saturated intrinsic permeabi 1 i ty(cm2) 
c c»pore disconectedness index 
c Zr=surface layer tkickness(cra) 
c Mo*vegetation cover 
c Kv=plant coefficient 

c k=parameter of gamma distibuted storm depths 
e Lamda=parameter of gamma distributed storm depths 

C n**’¥******^****m^t****'>'^i‘*mm** ******** *-****m*****’^***i'**t*iit***m 

min,mo,m,n,nu,kl ,mtb,mtr,mh, in 
T<jal sjk(20) .yi(20).soj(20) ,a77(20) ,b77(20) ,b78(20) 

J>®a1 da(365) ,SKP( 363), St (365). b79(20) .379(20) ,ys(20),yg( 20) 
real a78(20).day(365) 

f ii(d,so)=1 ./(d*(i.-so)**( 1.45-.0375*d)+5./3.) 
external plot_$setup (descriptors) 
external p1ot_$scale (descr iptors) 
external plot_ (descriptors) 
character=i‘ 10 xaxis.yaxis 
f i (em)“10.*t( .66+.5,5/em+. 14/em**2. ) 
kl1 = l 
ran=1. 

print, 'To use Manabes parameterization type 2 , otherwise 1' 
input ,mnb 

if(mnb.eq.l) go to 3020 

print, 'Input the initial soil moisture so' 
input, so 
go to 3021 

3020 print, 'Input the average annual soil moisture so' 
input, so 

print, 'Input Time step (in days) ' 
input, tis 
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C NR=Number of rainstorm events you want to generate 

3021 print, 'Input NR' 
i nput , NR 

print, 'Input storm properties k and Lamda' 
input.xk.aml 

print, 'If you want a simulation period greater than mtau print 2, otherwise 1' 
input.SPE 

if(SPE.eq.2) go to 4036 
go to 4037 

4056 print, 'Input the number of simulation periods' 
input,NSPE 

4057 if(mnb.eq,2) go to 3040 

print, 'For daily fluxes type l.for half hour fluxes type 2 ' 
input, fl 

print, 'To plot S(t) type 2, otherwise 1' 
input, lot 

print, 'For cumulative fluxes after each storm and interstorm period type 2, otherwise I' 
input,ucu 

print, 'To plot S(t) for different values of Zr type 2 .otherwise 1' 
input.szr 

print, 'To print the cumulative fluxes only at the end of the rainy season type 2 . otherwise 
input.fcu 

3040 print, 'To print the rainstorm events type 2 . otherwise 1' 

input.rae 

n«i 

3003 print, 'epr,mtb,mtr,mpa,mtaU,ta,mnu,n' 
input .epr.mtb, mtr.mpa, mtau. ta,mnu,n 
if(mnb.eq.2) go to 3022 

2020 if(SPE.eq.l) go to 2021 

if (SPE,eq.2.and.ran.eq.2) go to 4053 

2021 print, 'Mo, Kv,k1,c,Zr' 
input. vg,vk,k1 ,cs, zr 
Iffvg.eq. 1) stop 
if(ran.eq.2) go to 3004 
if(dif.eq.2) go to 3004 

C d(s)=evapotranspiration efficiency function 
C Ys(s)=surf ace runoff function 
C Yg(s)=ground water percollation function 

1000 print, 'To plot iJ(s) and y(s) type 2 , otherwise 1' 
input.pl - 

if(pl.eq.l) go to 3004 
1f(kl1.eq.2) go to 3004 

print, 'To draw different curves for J(s) for different climates type 2. otherwise 1' 

Input, dif 

double precision sum 1 .meant ,mean2.mean3,B28 
double precision sum2 
double precision sum3 

3004 if(ran. eq.2) go to 807 

3022 if(rae.eq.i) go to 42 

print.'STORM DEPTH STORM DURATION TIME BETWEEN ' 
print,' (cm) (days) (days) ' 

42 i 1-1 

C GENERATION OF RAI STORM EVENTS 
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C R1(I)=storm depth(cm) 

C R2(I)=storm duration(days) 

C R3( I )= 1 nterstorm duratlon(days) 

real R ( 3000 ) . WK ( 6000 ) , R 1 ( 3000 ) . R2 ( 3000 ) . R3 ( 3000 ) 
dcsuble preeicion DSEEO 
OSEED“ 123765. ODO 
A=xk 
B- 1 . /ami 

call ggamr(DSEED,A,NR,Wk.R) 
do S 1=1 ,NR 
B(I)=B*R(I) 

5 continue 
do 41 I = 1 . NR 
R1(I)=R(I) 

4 1 cont i nue 

0SEED=3478758.0D0 

A»1. 

B*mtr 

call ggamr(DSEED,A,NR,WK.R) 
do 7 1*1 , NR 
R(I)=B*R(I) 

7 continue 
do 21 1=1, NR 
R2(I)*R(I) 

21 continue 
OSEED=649353.0D0 
A* 1 . 

B*mtb 

call ggamr(OSEED,A,NR,WK,R) 
do 9 I = 1 , NR 
R(I)=B*R(I) 

9 continue 
do 30 1=1 .NR 
R3(I)=R(I) 

30 continue 
if(ran.eq.2) go. to 807 
1f(rae.eq.1) go to 3023 
go to 3024 

3023 If (mnb.eq.2) go to 3025 
go to 807 

3024 do 1 1 1*1 , NR 
wr1te(5,17) R1(I).R2(I),R3(1) 

17 format(f 10.6, 4x,f 10.6. 4x.f 10.6) 

1 1 cont 1 nue 

807 m=2./(cs-3. ) 
d*cs-1 ./m-1 
dE=2. l./m 
f1ed=f1e(dE) 

C if >r.iM4i*m**5|t*.t*** + + *#*>t»*5*r*s*i**%%5*i*5(‘1f* 3 *3t«;*j)e«**^3ft***^* 

C COMPUTE WATER CONSTANTS 
ca 1 1 WATCN ( ta , sut , nu , gamsw ) 

Q ifc * 4c ic 4c it ic 4: i^iciciciticiticiciiiltititiri^ 

c COMPUTE CLIMATIC PARAMETERS 

Q 4ciMft4i4cit4iiMt(ir4t4t4fiiic4c4e4i4r4M|M|(4(4i4t4r4(4r4iifi[.ici(;>*ci(.>4i4t4it i(*ic 4rA4c4'4citi'4i*'<ic*ir*ic4cic4c9tc4c‘4;4r 
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del ta»1 ./mtr 

mh^rnpa/ ( mtau/( mtb+mtr ) ) 

amnu=mtau/ ( rntb+rntr ) 

m I “mh/mtr 

ota= 1 ./mh 

alpha^l ./mi 

pl=3. 14159 

beta= 1 ,/nttb 

C COMPUTE DERIVATIVE OF d WITH RESPECT TO so 

den»( 1,425-0. 037S*d) 

If(pl.eq.l) go to 802 

k*0 

so»0. 

805 so»so+0.05 
go to 802 

802 ds»( 1 . -so)**den 

dds>ds*d 

deno»dds+(5./3. ) 

denom«'deno**(4./3. ) 

soo*1.-so 

so 1 »soo* » ( -4 . /3 . ) 

denos=2*soo*deno 

dt»(2. 425-0. 0375*d) 

so2“seo**dt 

deos1=so2’*d*den 

nom=-denos-deos 1 

nom1’'nom*so1 

dep=nom1/(denom*3) 

f 1c*f 1(m) 

si 1*sqrt(n/(kl*f lc))’t>sut/gamsv/ 
si 11“sl 1»so**( ^1 ./m) 
bk 1 ■ k 1 *gamsw/r )0 

8lgc=n*eta^*2 . *bk1*s1 1/(p1 *m*del ta) *72000. 

s1gc1*s1gc* *0.3333333 

ders 1 g=s i gc 1 *der 

s1a=5*n*bk1*86400*s1 1/(3*m*p1 ) 

slgma®(sigc/deno*( 1 . -so)**2. >**.333333 

g*alpha*bk1 *86400* . 5*( 1 .+so**cs) 

gl=al og10( sigma) 

xp=( 1 .766*g1 ) + (0,9S0*(g1**2, )) 

xp1»- .806-xp 

CSI»10. **xp 1 
xp2=( 1 ,96*gl )+1 .766 
U*-derslg*xp2/sigma 

co*alpha*86400*bk1/2. *cs*so**(cs- 1 . ) 

co1=U-co 

C2=*co1 *CSI *exp( -g) 

C38»mtau*86400*bk 1*cs/mpa*so**(cs- 1 . ) 

C3=C33/2, 

If(vg.eq.O) go to PO 
go to 90 


80' E*2.*beta*n*bkl*s1 1*f Ied/(p1*m*epr**2. )*86400*so**(d+2. ) 
If (E.ge.88. ) E=33, . 

z1=( 1,+E*sqnt{2, ))*exp(-E) 

22=gamma( 1 .5)-gamt( 1 .5, E) 
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z2az2’*sqrt(2, *E) 

sj = 1 , -z1+z2 

f f (pi .eq. 1 ) go to 803 

k = k+ 1 

sjk(k)=sj 

if(k,eq.20) go to 804 
go to 805 

803 ag^gammaCl .5)-gamt( 1.5.E) 
g1=exp(-E)*sqrt(2.) 
g2«£*sqrt(2, )+l . 
g2=g2*exp( -E) 

g3*ag*sqrt(2. )/(2.*sqrt(E)) 

g4=*axp(-E)*sqrt(E)*sqrt(2*E) 

gg»-g1+g2+g3-g4 

E1 1*2 . *beta'*n*bk1 *81 1 *f led/(p1 •mtepr*>»2. ) *86400 
E12*(d+2. )*soi'*(d+1 . ) 
deri j=gg*E1 1*E12 
C1«derl j 

<f (Cl . Ie.0.01) C1=0.0 
go to 100 

90 B=( 1.-vg)/( i. + (vg*vk)) 

B*B+(Vk*vg*-*2. )/(2.=f( 1 . + (vg*vk) )**2 . ) 

C*1 ./(2. *(vg*vk)»*2. ) 

E1«2. ♦beta*n»bk1*s1 1*f led/(p.i *in*epr*'*2.) *86400 

E*2. *beta*n*bk1 *sl 1*f1ed/(pf *m*epr**2 . )*86400*so**(d+2 . ) 

o1=B*( (vg*vk)+1 ) 

ol«-o1+sqrt(8*2.) 

ol 1=B*E*Sqrt(2.*B) 

o1»o1-o1 1 

o1*o1*exp(-B*E ) 

o1«o1*E1*(d+2. ) 

o1«o1*(so**(d+1 . ) ) 

o2*-vg*vk*C 

o2*o2+sqrt(2*C) 

o2-o2-(C*sqrt(2*C)*E) 

C88=C*E 

(F(C88.ge.88) C88*88. 
o2=o2*exp( -C88)*E1*(d+2 . ) 
o2“o2*(so**(d+1 . )) 

CE-C*E 

BE=B*E 

a1*(vg*vk)+1 . 
a2«E*sqrt(2. *B) 
a3-a1+a2 

If (BE.ge.aa. ) BE=88. 

a3»a3*exp(-BE) 

a4=vg*vk 

a4*a4+(E*sqrt(2.*C) ) 

1f(CE.ge.88. ) CE=88. ’ 
a4*a4*exp(-CE) 

a5*gamt(1 .5,CE)-gamt(1 .5,BE) 
aS*a5*sqrt ( 2 . *E) 
a6=a3-a4-a5 

a6=aS*( 1 . -vg)/( 1 . ~vg+(vg*vk) ) 
sj*1 . -a6 

if(pl.eq. 1) go to 806 

k»k+1 

sjk(k)*sj 

1f(k.eq.20) go to 804 
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go to 805 

806 o3=gamt(l.S.CE)-gamt(l.5,BE) 
o3“o3*?,qrt(2. *E1 ) 

03 -( I ,+cl/2. )*o3» (so**(d/2. ) ) 
o31»'C«El »(so-«(d+2. )) 
o31»(C*«1 .5)*exp(o31 ) 
o32»-B’i'E1»(so«*(d+2. ) ) 
o32*(3** 1 .5)*exp(o32) 
o33“o31 -o32 
033»033*(Ei*«1.5) 
o33-o33*(2.+d) 
o33»o33*(so**( ( 1 .5’»d)+2. ) ) 
o33=o33»sqrt ( 2 . *E } 
o3*o3+o33 
derj =o1 -o2-o3 
dopj=darj*( I. -vg) 
derj *'derj A (vg*vk)+1 . -vg) 

C 1 »der j 

<f(C1.1e.0.01) C1=0.0 
B2Q’'mtau»bk 1 tSo400/nipa*so*-''CS 

C *♦»♦*)*«*•** ♦•***0****l(l****1t*******»***#«»***««**** 

C Cl •'Derivative of J with respect to s 
C C2’-Derivati ve of Ys with respect to s 
C C38=»Derivativa of Yg with respect to s 
C sjoJCso) 

C sill^psi evaluated at so 

C bM^saturated hydraulic conductivity (cm/sec) 

C i **•'»•«*****» *Sk«*««*«*1* m 

100 prfnt10i-C1.C2,C38,sJ,sl11.bk1 

101 format(3hCT».f 10.6, 4x.3hC2«.f 10.6, 4x,3hC3»,f 10,6. 4x.2hU=. f 10.6, 4x.3hMH».f 10.2, 4x.f 20. 
SK“SO 

804 p*mpa/(mnu*mtr) 

Cl=C1*epr 

{f(sj .ge.0,99) sJ-1,0 
B1*sJ *epr 

If(pl.eq.l) go to 808 

so»0. 

k»0 

811 SO=SO+0.05 
ds»d*( 1 , -so)**den 
denoisds+Cs./S. ) 

sigma=(sigc/deno*(i . -so)**2. )**. 333333 
808 B22=s1gfJia*’»'( -sigma) 
sigm=slgma+1 . 

B22=B?'**gamma ( s 1 gm ) 

B2='.r^ii. exp(-g-(2*sigma)) 

B28-ntau*bk 1 '•86400/mp3*so’*'t>cs 
B4'=B2*p 

B5»'B2a*p’i'mnu*mtr/nitau 

if(pl.eq.l) go to 809 

k=k+1 

ys(k)=84 

yg(k)=B5 

soj(k)*'so 

if(k.eq.20) go to 810 
go to 8 1 1 
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809 1f{ucu.eq.2) go to 1816 

print,' S(t) i(cm/day) Et(cm/day) y ield(ctn/day) 

go to 1815 

C CALCULATE THE SOIL MOISTURE CONCENTRATION AND THE 
C CUMULATIVE EVAPORATION AND YIELD AT THE END OF 
C EVERY RAISTORM AND INTERSTORM PERIOD 

1816 print. 'SOIL. MOIST, CUM.EVAP, CUM. YIELD' 

1815 if(pl.eq.l) go to 812 

810 1f(k11.eq,2) go to 3001 

Q **»*ik-4 M* *‘1' *****’***<***»'********* ****** 

C PLOT d VERSUS s 

Q m*m*** ****************************************** *************** 

call plot_$setup( ' 's' , 'J' , 1 ,0,0.0) 

call plot $scate(0. , 1 . ,0. , 1 . ) 

3001 1=0 
kl l»2 

do 813 j=1,20 
1*1 + 1 

b77(l)=sjk(j ) 
a77( 1)=soj(j ) 

813 ccntinuQ 

call plot_(a77,b77,20, 1. ' ') 

If(dlf.eq.l) go to 3002 
go to 3003 

3002 read(5, ) 

C ************************ *************************************** 

C PLOT Ys AND Yg VERSUS s 

q *********** **************************************** ********** 

call plot_$setup( ' '.'SOIL MOISTURE 'SURFACE RUNOFF' . 1 ,0,0.0) 
call plot SscaletO. . 1. ,0: .2. ) 

1-0 

do 814 J=1,20 
1-1 + 1 

b78( 1 )=ys< j ) 
a78(1)*soj(j) 

814 continue 

call p1ot_(a73,b78,20, 1 , ' ') 
rcad(5, ) 

call plot_$setup.(' ','SOIL MOISTURE 'GROUNDWATER RUNOFF 1 .0,0.0) 

call plot $scale(0. , 1 . ,0. .2. ) 

1-0 

do 834 j=1 ,20 
1-1 + 1 

b79(i)=yg(J) 
a79( 1 )=soj ( j ) 

834 continue ■ 

call plot_(a79.b79,20, 1 , ' ') 
go to 1000 

812 1f(szr.eq. 1) go to 817 

do 2001 n- 1,2 

print, 'Input Zr(cm)' 

Input, zr 
817 a»n*zr 
Dt»t1s 


DAY ' 
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LS»0 

4053 If (SPE.eq.2,and.ran.eq.2) go to 4054 
1*0 

go to 4055 

4054 LS=LS+1 
If(LS.ge.NSPE) stop 

4055 LM=0 

K»0 

KP=0 

tsc^O-O 

SK3-0.0 

SK2»0.0 

LMM»0 

y f aldceO.O 

evapc^O.O 

400 if(ucu.eq.l) go to 401 
if(szr.eq.2) go to 401 
If (fCu.Bq.2) go to 40.1 
wr1te(6, 1701 )5K,avapc,y1e1dc 
1701 fdrmat(f8.S.4x,f3.5,4x,f8.5) 

C CALCULATE THE VALUE OF SOIL MOISTURE EVERY HALF HOUR 
C DURING A PRECIPITATION EVENT 


401 Dt1=0. 

yt-O.O ' 

slal^sla^f i i (d.SK) 
s1a2»2*( 1 , -SK)*sqrt(sla1) 

Ao«bk1+Se40\,‘'2. 
if(SK.le.O) go to 215 
ao1“Ao*( 1 .+<SK**cs)) 
go to 216 

215 ao1»Ao 

216 i«It1 
r2=R2(I) 
in=R1(I)/r2 
To1=»2»1it»( In-aol ) 
to2=sia2**2./To1 
to3«2. •( in-aol ) 
to4=1 .+(ao1/to3) 

To=to2*to4 

300 Dt1=Dt1+Dt 

if (Dtl .ge.r2) go to 200 

LM=LM+1 

If (Dt1 .ge.To) yt=1 
if (SK2. 1t.SK3) yt=0.0 

SK1=SK+( in-p'*((B2’*'yt)+{B28=»mnu*mtr/mtau) )-p*(SK-so)»( (C2*yt )+(C3*milu*mtr/mtau) ) )*Dt/a 

SK2=SK1 

SK3=SK 

If (SKI. ge. 0.999) go to 211 
go to 212 

211 SK1=0.999 
yield® In 

y1e1dc®yieldc+( in*t1s) 
go to 213 

212 y ield»p*( (B2’*yt)+(PT8*mnu*mtr/mtau) )+p*(SK-so)*( (C2*yt)+(C3*mnui'mtr/mtau) ) 
yfe1dc»yie1dc+(yield*t1s) 
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213 SK»SK1 

1 f (f 1 .eq. 1 ) go to 250 
<f(srr.eq.2) go to 300 
wr1te(6,210) SX. In, yield 
210 f ontat( f 3 . 5 , 4x . f 8 . 5 , 22x. f 8 . 5) 
qo to 300 

250 tis5=1,/tis 

If (LM.ge. tiss) go to 251 
go to 300 

251 LM«0 
LMM»LMM+1 

If (UMM. It.ratau) go to 907 

wri te(6,903) SK.evapc.yieldc.I ,tsc 

908 forniat(f8.5.4x.fa.5,4x,f8.S,4x, 1 10.4x,f 10.3) 

go to 900 

907 KPaKP+1 

SKP(KP)=SK 

da(KP)*LMM 

1f(ucu.eq.2) go to 300 
1f(szr.eq.2) go to 300 
wrfte(6.252) SK. In. yield, LMM 

252 format(f8.5.4x.f8.5,22x,f8.5,9x. 15) 
go to 300 

200 if(ucu.eq.l) go to 201 
lf(fcu.eq.2) go to 201 

wr1te(6, 1700) SK.yleldc.yt 
1700 format(f8.5, 16x,f8.5,4x,f3. 1) 

Q ■^*nm*ifx*^4iM**<tt*‘***m**<ititi**!**mm*** **««^*tt>l<**«****ii«ii'^*>l<*****>t'f<>«li^ 

C CALCULATE THE VALUE OF SOIL MOISTURE EVERY HALF HOUR 
C DURING AN INTERSTORM PERIOD 

Q *^m^*»**rniit’Himit***m**m*****ii0**»***»m»**»‘*f*»m***********^m*x- 

201 Dt1=0, 

500 Dt1=0t1+Dt 
r3=R3(I) 

If (Oti .ge.r3) go to 400 
LM*LM+ 1 

evap^B 1+(C1*(SK-so) ) 

If(ev^p.ge.epr) go to 600 
tsc»tsc+t 1 s 
•vapp=evap/epr 
If ievapp. le.vg) go to 701 

SK1»SK-(evap+(B28*p*mnu*rai;3-/mtauiM(c3’»p-*'mnu*K t?‘*(SK'V5o)/ffltau) )i*Dt/a 

If(SKI.Ie.O) SK1=*0.0 

•vapc=evapc+ (evap* t 1 s ) 

go to 700 

600 evap=epr 

evapc=evapc+ ( evap* 1 1 s ) 

SK 1 “SK- ( epr*D t/ a ) - ( ( B28 *p*mnu*mtr/mtau ) + ( C3 *p*mnu*nitr* ( SK -so)/mtau ) ) *0t/a 

If(SKI.le.O) SK1=O.0 

go to 700 

701 evap=epr*vg 

6vapc“evapc+(evap»t1s) 

SKl“SK-(evap*Ot/a)-( (828*p*mnu*mtp/mtau)+(C3*p*mnu+mtr*(SK-so)/mtau) )*Dt/a 
If(SKl.le.O) SK1“0.0 

700 y ield“(B28*p*mnu*mtr/mtau)+(C3*p*ninu*mtr*(SK-so)/mtau) 

If (yield, le. 0.0000001) y lei 0*0.0000001 
y1eldc*yleldc+(y1e]d*tls) 

SK-SK1 

if(f 1 .eq. 1) go to 750 
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if ('szr.eq.a) go to 757 
writts(6,220) SK.evap, yield 
220 for.t^at(f8.5. I6x, f8 . S . lOx , f 8 . 5) 

757 K“K+1 

if (K.ge. ■5000) stop 
go to 500 

750 t1ss»1 ./tis 
if (LM,ge. tiss) go to 75 t 

go to 500 

751 LM=0 
LMM»LMM+ 1 

f f (LMM. le.ratau) go to 901 
wri te(S,905) SK.evapc.yieldc. I , tsc 
905 forfnat(f8,5,4x,fa,5,4x,f8.5,4x. i 10,4x,f 10.3) 
go to 900 
901 KP»KP+1 
SKP(KP)=SK 
da(KP)=LMM 

1f(ucu.eq.2) go to 500 
if(szr.eq.2) go to 500 
write(6.752) SK.evap. yield, LMM 

752 format(f8.5, ISx.fS.S, t0X,f8.5,9x, i5) 
go to 500 

900 if(szr.eq.2) go to 2003 
if(ran.eq.2) go to 5000 

C *m*ii^*m’tm**^‘****'***********’tr**************»******t»**i*****’t0** 

C CALCULATE THE STATISTICAL PROPERTIES OF THE GENERATED 
C RAINSTORM EVENTS 

C * *ifim** ******* m**!* ********•**•'*’•*•*’•'*’****•*'* ********** 

print, 'Statistical properties of the simulated rainstorm characteristics' 

5000 5um1»0.0D0 
sum2=0.000 
sumSaO.ODO 

do 1001 IL*1,I 
sum1»sum1+R1 (IL) 
sgm2=sum2+R2(IL) 
sunj3=sum3+R3( IL) 

1001 continue 
mean1=sum1/(f ioat(I ) ) 
mean2 >suhn2 / (float(I)) 
me3n3=sum3/(f loat(I ) ) 
var1’=0.0 

var2“0.0 
var3*0.0 
do 1002 IL=1,I 

var1=var1 + ((R1(IL)-mean1)*»2. ) 
var2=var2+( (R2{ rL)-mean2)*f2. ) 
van3=var3+( (R3(IL) -mean3)**2, ) 

1002 continue 
variiavar 1/f loat(I- 1) 
varl2=var2/f loat(I-1 ) 
vari3avar3/f 1 oat(I- 1 ) 

If(ran.eq.l) go to 5001 
NSPP-NSPE-1 

if (ran.eq.2.and.LS.eq.NSPP) go to 5001 
go to 2020 

5001 print. 'AVER. h(cm) AVER. tr(days) AVER. tb(days) ' 

wrlte(6, 1003) meani ,mean2,mean3 

1003 formac(f10.6.6x.f 10.6,6x.f 10,6) 


C - 
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print,' VAR.h VAft.tr VAR.tb' 

print tOOA.vari 1 .vari2,variti 
1004 format(f8.2.4x,f8.2. 10X,f3.2) 
ran=2 , 

■2031 if(lot.eq.2) go to 2030 
go to 2020 
2030 read(5,) 

2003 If(n.gt.l) go to 2003 

C PLOT THE SOIL MOISTURE CONCENTRATION WITHIN THE 
C LAVER OF THICKNESS Zr VERSUS TIME DURING THE 
C RAINY SEASON LENGTH 

C If «»«%**«** «*» *4i#i«lk*«i««**«**4'«-ltt'V****’«:***4<*** 

call plot_$setup( ' 'DAYS' . 'SOIL MOISTURE' , 1 ,0.0.0) 
cal 1 plot_?scals( 1 , .220. .6. , 1 . ) 

2003 1=0 

do 910 J-1.LMM 
1 = 1 + 1 

St(l)=SKP(j) 
day( 1 )=da( j ) 

910 continue 

If(ll.eq.l) go to 2004 

lf( n .eq.2) go to 2005 

2005 call plot__{day,st,mtau.3, ' i ' ) 

go to 2001 

2004 call plot_(day,5t,mtau, 1 , ' ') 
ir(szr.eq.l) go to 2000 

2001 continue 

C CALCULATE THE MOISTURE FLUXES USING MANABE'S PARAMETERIZATION 
C 4i**4i*<<t«4t4i**;^4t4r***«*j|ti«i*4r#*4Mtt4i««4!***«****»*4>**«*4<ii«****»4i«*i«**»«* 

3025 1f(mnb.eq. 1) go to 2000 

print. 'S(t) CUM.EVAP. CUM. YIELD' 

SK*so 

0t«1./48. 

I«0 

yleldc=0.0 

evapc=0.0 

Dt11=0.0 

3031 wrlte(6,3033) SK.evapc.yletdc 
3033 format(f8.5.4x.f8.5.4x,f8.S) 

Dt1=0.0 
I«I + 1 
r2=ri2(l) 
ln»Rl(l)/r2 
30;i:8 Dt1=0t1+Dt 
Dt:i 1=Dt1 1+Dt 

lf(SK.ge.0.42) go to 3029 
SKl»SK+ln*Dt/(n*100) 

SK=SK1 

If (Oti -ge.r2) go to 3027 
go to 3028 

3029 y1eld=( ln-epr)*Dt 
y1eldc*yleldc+yleld 
lf(Dt1 .ge.r2) go to 3027 
go to 3028 

3027 wrlte(6,3030) SK.evapc, yields 
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3030 f ormat (f 3 . 5 . 4x . f 8 . 5 , 4x . f 8 . 5 ) 

0t1“0.0 

r3=R3(I) 

3032 0t1-Dt1+0t 
Dt 1 I*0t1 1+Dt 
evap=epr 

Jf (SK. 1 1. 0.315) evap=»epr*SK/0.315 
SKlBSK-e'/ap^0t/(n*100) 
evapc*evapc+ ( evap*D t ) 

SK=«':,K1 

If (Dtll .ga.mtau) stop 

1 f (Dtl .ge, r3) go to 3031 
go to 3032 
2000 read (5.) 

stop 

end 

Q t^»lr^ml|^^^tl**a*^I^.^n*m*!*******'^**^^*ltt****■»*■****||^i*^:>>*m■•iL^i1lMMm********** 


subroutine WATCN( ta.sut.nu.gnmsw) 
real nu.nut 

dimension sutt( 1 i ) ,nut( 1 1 ) ,gamst( 1 1 ) 

data sutt/75. 6. 74.9, 74. 2. 73. 5. 72.80. 72. 1 .71.4.70.7.70.0.69.3.68.6/ 
data nut/17.93e-3, IS. 18e-3. 13.09e-3. 1 1 .44e-3. 10.03e-3.8.94e-3. 

& 8,e-3.7,2e-3.6.53e-3.5.97e-3,5.94e-3/ 

data gams t /0 .99987.0. 99999999 . 0 . 99973 .0. 999 13.0. 99823 . 0 . 99708 , 

8 0.99568 . O . 99406 .0 . 99225 , 0 . 99025 .0.98807/ 

If (ta.gt.50. )go to 10 
1ta=if ix(ta*.2)+1 
f rac»ta-f 1 oat(5*( i ta- 1 ) ) 

1ta1»ita+'< 

sut=(sutt( i ta1 )-sutt( 1 ta) )*0.2’*'f rac+sutt{ 1 ta) 
nu»(nut( ital )-nut( 1ta))*0.2*frac+nut( ita) 
gamsw=( (gamst( i tal )-gamst( ita) )*.2*f rac+gamstC 1 ta) )*980. 
return 

10 sut=sutt(11) 

nu«nut(11) 

gamswsganist( 1 1 ) 

return 

end 

c this function computes the gamma incomplete function 

function gamt(a.x) 

1f(x.eq.0)gp to 40 
if (x.gt. 100)go to 50 
sum=1 ./a 
an" 1.0 
old“sum 

33 o1d=ol d*x/(a+an) 
if (old/ sum-1 .6-6)20.10. 10 
10 an=an+1. 

Eum=5um+o1 d 
if(an-300. )33.33. 12 
12 continue 

20 xxx=(a-*alog(x)+a1og(sum)-x) 
if (xxx. 1 t . -80. )go to 40 
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9 afnt = lexp(xxx) ) 
go to 60 
40 gamt=0.0 
go to 60 
50 gamt=gamma(a) 

60 return 
end 

c This function computes the gamma function by a Stirling approx. 

function gamma(y) 
x=y+1 . 
pl«3. 14159 
St1r1=1 ./( 12.*x) 

Stir2=1 ./(288.*x**2.) 

Stir3=-139,/(5ia40. *x*+3. ) 

Stir4='-571 ./(2488320.1'X**4. ) 

St lr“1+stir1+st 1 r2+st1r3+stir4 
gamma=exp(-x)*x*t‘(x> .5)'*sqrt(2. *pi )*st ir/y 
end 

function f1e(d) 
dimension y(6) 

data y/0. 13,0.11 ,0.077,0-056,0.044,0.034/ 

If (d.gt.7 . )go to 10 
X»d- t . 

1«1f ix(x) 

fpac»x-f loat( 1 ) 

y1«alog(y( 1 )) 

y2*alog(y( i+1 ) ) 

f 1e=exp( (y2-y 1 )*frac+y 1 ) 

return 

10 fie-. 034 

return 

end 
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C THIS PROGRAM CALCULATES THE LATENT HEAT 
C FLUX OBTAINED FROM THE ANALYTICAL MODEl, 

C USING EXACTLY THE SAME PERIODIC ATMOSPHERIC 
C FORCING SPECIFIED BY MILLY AND EAGLESONi 1932 . 

C TR.279) FOR THE CLIMATE OF WINSLOW. ARIZONA . 

C THE FUNCTIONS OF HYDRAULIC AND VAPOR CONDUCTIVITY 
C USED ARE APPROXIMATIONS OF THOSE SPECIFIED 
C IN TR.279 FOR SILTY LOAM AND SAND, 

C THE MODEL WAS ALSO MODIFIED TO TAKE INTO ACCOUNT 
C VAPOR FLOW FOR THE SILTY LOAM. 

C **iMr«* + i*>**(<*****^.t. •1i'*’l' + **»****ll(*****,***«*»*»»*»4 *»»*** 

c 

C CLIMATIC AND SOIL PARAMETERS 
C 

C Gp=average potential evaporation rate( cm/day) 

C mtb=mean time between starms(days) 

C mtr-mean storm durat ion(days) 

C mt=mean rainy season length(days) 

C ta^average annual temperature(deg.Celcius) 

C mv=mean number of storms per year 
C mpa=mean annual precipitat ion(cm) 

C n-effectlve porosity of soil 
C k1*saturated intrinsic permeabi 1 i ty(cm2 ) 

C c»pore disconnectedness index 
C so=average annual soil moisture 
C yg=«average annual percolation rate(cm/day) 

C ys=average annual surface runnoff rate(cm/day) 

C mi=niean storm intensity (cm/day) 

C cl ,c2,c3=1 inear ization coefficients of the annual 
C water balance as obtained from the prog- 

C ram Taylor, fortran. 

C et=*actua1 average annual evapotranspiration rate(cm/day) 
C ad,aw=coef f icients of the albedo function as specified 
C in TR.279 

C Si='initial soil moisture 

C T1i=initial surface temperature(deg.celcius) 

C T2i=initia1 deep soil temperature(deg.celcius) 

C Zr=near surface storage depth(cm) 

C 

print , "ep.mtb.mtr ,mt . ta.mv.mpa” 

input .epr.mtb.mtr, mt , ta.mv.mpa 

pr int . "n. k 1 ,c . so, yg. ys ,mi . c 1 , c2 ,c3 . et ,ad. aw" 

input . un, aki , c , so . yg, ys, ami , cl . c2 . c3 , et , ad. aw 

100 print, "Si .TH ,T2i " 

input . sk , t Ik , t2k 

sini=sk 

print." Input Storage Depth 2r " 
input, Zr 

print, "If silty loam type 1,1f sand type 2" 

Input, soil 

if (soi 1 .eq. 1., ) spr=o. 0064426 

if (sol 1 .eq.2. ) spr =0.02438 19 

to=0.0 

prec=0.0 

toc=0.0 

tc=0.0 

t=o.o 

,k=0 

c3=c3/2. 

0t=0. 25/24. 
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m’‘0 

13 kt1=-74+(k*96) 
kt2=20+(k*96) 

( f ( tc. gt . k t1 .and. tc . 1 e . k t2) go to 11 
kt3=22+(k*96) 

15 1 f ( tc. gt .kt2.and. tc. Ie.kt3) go to 14 
k=k+1 

to=spr*2.+(sk**2. 165) 

to*86400*( to**2. ) 

to=to/(2. ♦{epr«*2. )) 

to=tO'»24 . 

pr Int , "To” 

wrlte(6,400) to 

400 format(f 10.4) 

toc=0.0 

prec= 1 . 

t1m=to 

go to 13 

C FINITE DIFFERENCE EQUATION FOR EVALUATING SOIL MOISTURE 
C DURING A DRY PERIOD. 

11 1*1 
yt*0.0 
cv=0.45 
toc*toc+0. 25 
1 f ( to.gt . toe) sic*1. 
if (to. le, toe) a1c=0. 
etaset+Ccl+epr^Csk-so) ) 

If (eta. 1e. 0.0) eta*0,0 
yaa=yg+( (sk-so)*c3*arni ) 
if (yga. le.0.0) yga=0.0 

if (to.ge. toe) eta=0.00 1 19*0.622*ch*Ua*3600. t-Ces-g .67)»24 ./ 1013 .25 
1 f ( to. 1 e. too. and.prec. eq. 1 . ) go to 600 
go to 601 

600 If (so1 1 .eq.2) go to 601 
th1*un*sk 

if (thi . le.O.O) go to 601 
y1m=(-30.*th1 )-5.5 
Os»(10.*»(-5.5)) 

De®De-( 10. **y1m) 

Oe=De/1.5 

De=De>^3600. + 1 . 85* ( th 1 ** ( - 1 .83 ) ) 

Se»2 . *th1* ( sqrt (De/3 . 14 ) ) 
t ini=tim+0. 25 
fet=Se/(2. •sqrt(tim)) 
fet=fet-(2.852e-9) 
fet=24.*fet 

1 f ( f et . 1 t . eta) pree*0.0 
i f ( f et . ge, eta . and. sk.ge.s1ni) eta=fet 

601 sk 1=sk+( (Dt/(un*Zr ) ) *( -eta-yga) ) 
if (ski . 1 t .0.0) sk1=0.0 
eta=eta*597./24. 

if (ski .eq. 0.0) eta=0.0 
go to 12 

C FINITE DIFFERENCE EQUATION FOR EVALUATING SOIL MOISTURE 
e DURING A RAINY PERIOD. 

14 1=0 
yt= 1 .0 

HP=0. 193*0. 99*Ta 
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CV= 1 . 

sk 1 =3k + ( (D t/(ut>'Zr ) ) *(4 .S32-yg-ys- ( (sk-so) •■amf * (c3+c2 ) ) ) ) 

sk 1 =sk 1 - ( Ep’Ot/Curr^Zr-oS? . ) ) 

ifCskI.gt. f .0) skl=1.0 

C Air? TEMPERATURE 

12 argst- 15 . 

arg=argT3. 14/ 12. 

Ta=25.8+(7 .8*cos(arg) ) 
c WIND SPEED 
argi =t - 18.5 
arg1=arg1 *3.14/12. 

Ua=3S0.+( 180. *cos(arg1 ) ) 

R 1 *392400. /(Ua* *2. ) 

R1*R1 *(Ta-t1k)/(ta+546.+t1k) 

1 f(R i . 1 t .0,0. and. R1 .ge. -0.0014) ch=0.0028 
if(Rl , le. -0.0014, and, Ri.ge. -0.0054) ch*0.0029 
lf(Ri . 1e. -0,0054, and, Ri .ge. -0.0105) ch*0.0030 
if (R1 . 1e. -0. 0105. and. RI ,ge. -0.0205) ch=0.0032 
lf(Ri , le. -0.0205. and. Ri ,ge. -0.0402) ch=0.0035 
if (RI . le. -0.0402. and.Ri .ge. -0.0793) ch»0.0039 
if (Ri . le. -0.0793.and.Ri .ge.-O. 1569) ch=0.0044 
if(Ri . le. -0. 1569.and.Ri .ge, -0.3119) ch*0.0052 
if(Ri. le. -0.3119) ch»0.0058 
if (Ri .ge.0.2) ch*0.0 

if (Ri .ge.0.0. and.Ri . 1 1.0.2) ch*0.00277>>'( ( 1 . -(R1/0.2))**2. ) 

C EVAPORATION RATE 

es=6. 1 1+(0.6102+t1k) 

Ep*0.001 19*0.522*ch*Ua«3600.*(es-9.67)/1013,25 

Ep*Ep*597. 

if(l .eq.O) EL=Ep 

if(eta.lt.Ep) th1=0.0 

i f (eta.ge.Ep) th1=un*sk 

if (1. eq.O) th1*Un^sk 

if(sic.eq.l) th1=un*sk 

if (thi .eq.O.O) eps*0.9 

if (thi igt.0.0) aps*0.95 

if (eta.ge.Ep. and. 1 ,eq. 1) EL=Ep 

if (eta- 1 t. Ep.and. 1 .eq. 1 ) EL=eta 

if (sic.eq. 1 . ) EL=Ep 

c ALBEDO 

thi 1=2. *th1 

if (th1 1 .gt.un) A=aw 

i f ( thi 1 . 1 e, un) A*ad+( (aw-ad) *th1 1/un) 

C SHORT WAVE RADIATION 

agu*t-12. 
agu=agu*3 , 14/ 12. 

sna=(s 1 n(0, 375) -*5 in(0. 61 1 ) ) + (cos(0. 375 ) *cos(0. 61 1 )*cos(agu) ) 

if (sna. le.0.0) go to 210 

a1=0. 128-(0.054*(alog10( 1 ./sna)) ) 

a 1n=a1 *^2 . 5/ sna 

if(a1n.ge.88. ) aln=88. 

if (aln. 1 e. -88 . ) a1n=-88. 

ae=exp(-a1n) 

vc*1 . -(0- 65*(cv**2 . ) ) 

wb=120. *sna*ae*vc 
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if (sna.gt.0,0) si'swb 
210 Iffsna le.O.O) si=»0.0 

C DOWN LONGWAVE RADIATION 
dl 1»{9. 37e-S)-(0.a26e- t0)*60, 
dl i=dl l*((Ta+273, )’' '■6. ) 
cvv=1 ,+(0. 17«cv**2. ) 
dl 1 =dH *cw 

C BACK LONG-WAVE RADIATION 
•Jl i=>(tlk+273. )»*4, 
ul i =>U1 I *(0.826e- 10) •60. *eps 
c SENSIBLE HEAT TRANSFER 
H=-0.001 19-*-0,2399*chi-Ua*3600. «(Ta-t1k) 

C FORCE-RESTORE METHOD FOR ESTIMATING 
C SURFACE TEMPERATURE, 
tlk=t1k+273. 
t2k=t2k+273. 

G=-ul i-EL'H>(HP*yt)+(( 1, -A)*sl )+(eps*dl 1 ) 

G«G-(0,99*Sic*( Mk-273)*Ep/"^97. ) 

rs=( ( 1 , -un) *2 .0e6) + {un*sin1 * 4 .2e6)+( (un-(un*sini))«1. 25e3) 
rs*rs/4,2e6 

If (sol 1 .eq. 1 ) aml = 1,3e-3 
if (sol 1 . eq. 2) aml=4.e-3 
ssk=aml/rs 
d2=20, •ssk*86400. 
d2=sqrt(d2) 

t2k»t2k+(0.25*G/(rs’*d2)) 

pri*( ( 1 . -un)*2.e6)+( th1+4.2e6)+( (un-thl )*1 , 25e3) 
pr1=prl/4.2s6 
ht=th1 

if (sol 1 .eq.2) go to 200 
If (ht.ge.0,4) alm=»3,75e-3 
if (ht. le.0.4.and.ht.ge.0. 3) alm=3.5e-3 
if (ht. 1e.0.3.and.ht.ge.O.2) alm=3.0e-3 
if (ht . le.0.2 .and.ht .ge,0, 1 ) alm*2.65e~3 
if (ht . le. O. 1 . and.ht ,ge. 0.05) alm= 1 . 5e-3 
if (ht. le.0.05) alnt=0.5e-3 
go to 300 

200 if (ht.ge.0.3) a1m=8.e-3 
if (ht . le. 0.3. and.ht .ge.0.2) a1m=7. 1e-3 
if (ht. 1e.0.2.and.ht.ge.0. 1) alm=6.4e-3 
if (ht. le.O. 1. and. ht.ge. 0.05) alm=5.e-3 
if (ht. le.0.05) a1m=2.e-3 
300 c1am=0.3*sqrt(alm*pM ) 
clam=clam+(0.7*sqrt(aml •rs)) 
ci1=2. *(sqrt(3'. 14/86400) ) /cl am 
C22=2. +3. 14/24. 

t1k=t 1k+(0.25*Cl l*G)-(0.25*C22*(t1k-t2k)) 
t1k=t1k-273. 
t2k=t2k-273. 

sk=sk1 
EL=EL+24. 

write(6,120) EL , t Ik , 1 . t2k .G.sk ,Ri 
120 format(f 10.4, 4x.f 10.4, 4x. 1 1 ,4x,f 10.4,4x,f 10.4,4x,f l0.4.4x,f 10.6) 
tC«tC+0.25 
t»t+0.25 

1f(t.eq 24.) t*0.0 
1 f O . eq . 1 ) go to 13 
go to 15 
end 
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